Bulls Eye for Athletes

Statistical analyses

Author
Affiliation

Gustaf Reinebo

Karolinska Institutet

Published

2025-02-13

1 Data preparation

The following data preparations are the basis for all analyses related to the Bulls Eye for Athletes (BEA) study. Code for scale step alterations and some of the modified data frames are presented next to their related sub-analysis to make the workflow more comprehensible.

Code
#loading bulls eye for athletes dataset
BullsEye_merged_data_final_without_comp_level <- read_excel("BullsEye_merged_data_final.xlsx")

rename <- dplyr::rename

#renaming bulls eye items from Swedish to English
BullsEye_merged_data_final_without_comp_level <-
  BullsEye_merged_data_final_without_comp_level %>%
  rename(
    BE_Competition = BE_Tavling,
    BE_Training = BE_Traning,
    BE_PreparationRecovery = BE_Forberedelser,
    BE_Life = BE_Livet,
    BE_Obstacles_Sport = BE_Hinder_Idrotten,
    BE_Obstacles_Life = BE_Hinder_Livet,
    BE_Competition_test2 = BE_Tavling_test2,
    BE_Training_test2 = BE_Traning_test2,
    BE_PreparationRecovery_test2 = BE_Forberedelser_test2,
    BE_Life_test2 = BE_Livet_test2,
    BE_Obstacles_Sport_test2 = BE_Hinder_Idrotten_test2,
    BE_Obstacles_Life_test2 = BE_Hinder_Livet_test2
  )

#adding competition level
comp_level_BE <- read_excel("comp_level_bulls eye.xlsx")

BullsEye_merged_data_final <- full_join(BullsEye_merged_data_final_without_comp_level,
                                        comp_level_BE,
                                        by = "Användarnamn")

#data corrections for BullsEye_merged_data_final
#age reported as year of birth
BullsEye_merged_data_final$Demografi.Alder[BullsEye_merged_data_final$Användarnamn == "1054fmsc" &
                                             BullsEye_merged_data_final$Demografi.Alder == 95] <- 24
BullsEye_merged_data_final$Demografi.Alder[BullsEye_merged_data_final$Användarnamn == "1180smpq" &
                                             BullsEye_merged_data_final$Demografi.Alder == 91] <- 28
BullsEye_merged_data_final$Demografi.Alder[BullsEye_merged_data_final$Användarnamn == "1183qhst" &
                                             BullsEye_merged_data_final$Demografi.Alder == 93] <- 26

#participant stopped practice sport by test2 so values at test2 therefore not included
variables_to_na_1086fcfd <- c(
  "BE_Competition_test2",
  "BE_Training_test2",
  "BE_PreparationRecovery_test2",
  "BE_Life_test2",
  "BE_Obstacles_Sport_test2",
  "BE_Obstacles_Life_test2",
  "SMS-2.1_test2",
  "SMS-2.2_test2",
  "SMS-2.3_test2",
  "SMS-2.4_test2",
  "SMS-2.5_test2",
  "SMS-2.6_test2",
  "SMS-2.7_test2",
  "SMS-2.8_test2",
  "SMS-2.9_test2",
  "SMS-2.10_test2",
  "SMS-2.11_test2",
  "SMS-2.12_test2",
  "SMS-2.13_test2",
  "SMS-2.14_test2",
  "SMS-2.15_test2",
  "SMS-2.16_test2",
  "SMS-2.17_test2",
  "SMS-2.18_test2"
)

BullsEye_merged_data_final[BullsEye_merged_data_final$Användarnamn == "1086fcfd", variables_to_na_1086fcfd] <- NA

#rename Valuing Questionnaire items abbrevation from VLQ to VQ
BullsEye_merged_data_final <- BullsEye_merged_data_final %>%
  rename_with( ~ gsub("^VLQ", "VQ", .), starts_with("VLQ"))


#create date and time difference variables, for survey 1 and survey 2
# Ensure the variables are in POSIXct format
BullsEye_merged_data_final$Färdig_surveys_test1 <- as.POSIXct(BullsEye_merged_data_final$Färdig_surveys_test1, format = "%Y-%m-%d %H:%M")
BullsEye_merged_data_final$Färdig_surveys_test2 <- as.POSIXct(BullsEye_merged_data_final$Färdig_surveys_test2, format = "%Y-%m-%d %H:%M")

# Calculate the time difference in hours (for survey 1 and survey 2)
BullsEye_merged_data_final$time_latency_survey_1_2_hours <- as.numeric(
  difftime(
    BullsEye_merged_data_final$Färdig_surveys_test2,
    BullsEye_merged_data_final$Färdig_surveys_test1,
    units = "hours"
  )
)

# Calculate the time difference in days (for survey 1 and survey 2)
BullsEye_merged_data_final$time_latency_survey_1_2_days <- as.numeric(
  difftime(
    BullsEye_merged_data_final$Färdig_surveys_test2,
    BullsEye_merged_data_final$Färdig_surveys_test1,
    units = "days"
  )
)



#create date and time difference variables, for survey 1 and bulls eye 1
# Ensure the variables are in POSIXct format
BullsEye_merged_data_final$Färdig_surveys_test1 <- as.POSIXct(BullsEye_merged_data_final$Färdig_surveys_test1, format = "%Y-%m-%d %H:%M")
BullsEye_merged_data_final$Skapat <- as.POSIXct(BullsEye_merged_data_final$Skapat, format = "%Y-%m-%d %H:%M:%S")

# Calculate the time difference in hours (for survey 1 and bulls eye 1)
BullsEye_merged_data_final$time_latency_survey1_bulls_eye1_hours <- as.numeric(
  difftime(
    BullsEye_merged_data_final$Skapat,
    BullsEye_merged_data_final$Färdig_surveys_test1,
    units = "hours"
  )
)

# Calculate the time difference in days (for survey 1 and bulls eye 1)
BullsEye_merged_data_final$time_latency_survey1_bulls_eye1_days <- as.numeric(
  difftime(
    BullsEye_merged_data_final$Skapat,
    BullsEye_merged_data_final$Färdig_surveys_test1,
    units = "days"
  )
)


#create date and time difference variables, for bulls eye 1 and bulls eye 2
# Ensure the variables are in POSIXct format
BullsEye_merged_data_final$Skapat <- as.POSIXct(BullsEye_merged_data_final$Skapat, format = "%Y-%m-%d %H:%M:%S")
BullsEye_merged_data_final$Skapat_test2 <- as.POSIXct(BullsEye_merged_data_final$Skapat_test2, format = "%Y-%m-%d %H:%M:%S")

# Calculate the time difference in hours (for bulls eye 1 and bulls eye 2)
BullsEye_merged_data_final$time_latency_bulls_eye1_2_hours <- as.numeric(
  difftime(
    BullsEye_merged_data_final$Skapat_test2,
    BullsEye_merged_data_final$Skapat,
    units = "hours"
  )
)

# Calculate the time difference in days (for bulls eye 1 and bulls eye 2)
BullsEye_merged_data_final$time_latency_bulls_eye1_2_days <- as.numeric(
  difftime(
    BullsEye_merged_data_final$Skapat_test2,
    BullsEye_merged_data_final$Skapat,
    units = "days"
  )
)


#have checked and all date and time variables are correct 

filter <- dplyr::filter

#hockeystudien bulls eye data:
sheet1 <- read_excel("hockeystudien intervention_bulls eye.xlsx", sheet = 1)
sheet1_filtered <- sheet1 %>%
  filter(`Instansens namn` == "Vecka 1") %>%
  distinct(`Användarnamn`, .keep_all = TRUE)

sheet2 <- read_excel("hockeystudien intervention_bulls eye.xlsx", sheet = 2)
sheet2_filtered <- sheet2 %>%
  filter(`Instansens namn` == "Vecka 1") %>%
  distinct(`Användarnamn`, .keep_all = TRUE)

merged_data_hockeystudien_BE <- full_join(sheet1_filtered, sheet2_filtered, by = "Användarnamn")

#hockeystudien demographics
hockeystudien_intervention_demographics <- read_excel("hockeystudien intervention_demographics.xlsx", sheet = 1)

#corecting age for a participant in the df
hockeystudien_intervention_demographics$Demografi.Alder[hockeystudien_intervention_demographics$Användarnamn == "1077dhrv" &
                                                          hockeystudien_intervention_demographics$Demografi.Alder == 6] <- 26


#merge df BE and demographics from hockeystudien
merged_data_hockeystudien_BE_demographics <- full_join(merged_data_hockeystudien_BE,
                                                       hockeystudien_intervention_demographics,
                                                       by = "Användarnamn")

#following in the hockeystudien study already exists in the bulls eye study, therefore mask:
#mask 1029spns, 1061bgmb, 1043gsnq, 1046pkbz, 1062xhqd, 1070mfqq

usernames_to_exclude <- c("1029spns",
                          "1061bgmb",
                          "1043gsnq",
                          "1046pkbz",
                          "1062xhqd",
                          "1070mfqq")

filtered_data <- merged_data_hockeystudien_BE_demographics %>%
  filter(!(`Användarnamn` %in% usernames_to_exclude))


hockeystudien_bulls_eye_modified_data <- filtered_data %>%
  mutate(`Användarnamn` = paste0(`Användarnamn`, "_hockeystudien")) %>%  # Append "_hockeystudien" to each value
  rename(
    BE_Competition = Tavling,
    BE_Training = Traning,
    BE_PreparationRecovery = Forberedelser,
    BE_Life = Livet,
    BE_Obstacles_Sport = HinderSkala1,
    BE_Obstacles_Life = HinderSkala2
  )




#preparations for merging datasets into superfinal version (bulls eye + hockeystudien_bulls eye)

hockeystudien_bulls_eye_modified_data <- hockeystudien_bulls_eye_modified_data %>%
  mutate(from_hockeystudien = 1)

BullsEye_merged_data_final <- BullsEye_merged_data_final %>%
  mutate(from_hockeystudien = 0)

bulls_eye_merged_superfinal_data <- bind_rows(hockeystudien_bulls_eye_modified_data,
                                              BullsEye_merged_data_final) #use for analyses

#BE items (dartboard), adjust (-1) to enable Rasch analysis
BE_adjust <- c(
  "BE_Competition",
  "BE_Training",
  "BE_PreparationRecovery",
  "BE_Life",
  "BE_Competition_test2",
  "BE_Training_test2",
  "BE_PreparationRecovery_test2",
  "BE_Life_test2"
)

bulls_eye_merged_superfinal_data[, BE_adjust] <- lapply(bulls_eye_merged_superfinal_data[, BE_adjust], function(x)
  x - 1)

#"Hinder", reverse items (low value = low in values-based living and compatible with other values items) and adjust (-1) to enable Rasch analysis
Hinder_reverse_and_adjust <- c(
  "BE_Obstacles_Sport",
  "BE_Obstacles_Life",
  "BE_Obstacles_Sport_test2",
  "BE_Obstacles_Life_test2"
)

bulls_eye_merged_superfinal_data[, Hinder_reverse_and_adjust] <- lapply(bulls_eye_merged_superfinal_data[, Hinder_reverse_and_adjust], function(x)
  8 - x - 1)


#ADAQ, reverse (higher value = higher PF) and adjust (-1) to enable Rasch analysis
ADAQ_reverse_and_adjust <- c("ADAQ.1",
                             "ADAQ.2",
                             "ADAQ.3",
                             "ADAQ.4",
                             "ADAQ.5",
                             "ADAQ.6",
                             "ADAQ.7")

bulls_eye_merged_superfinal_data[, ADAQ_reverse_and_adjust] <- lapply(bulls_eye_merged_superfinal_data[, ADAQ_reverse_and_adjust], function(x)
  8 - x - 1)


#SWLS, adjust (-1) to enable Rasch analysis
SWLS_adjust <- c("SWLS.1", "SWLS.2", "SWLS.3", "SWLS.4", "SWLS.5")

bulls_eye_merged_superfinal_data[, SWLS_adjust] <- lapply(bulls_eye_merged_superfinal_data[, SWLS_adjust], function(x)
  x - 1)

#VQ scale already starts at zero and does therefore not need to be adjusted for Rasch analysis

2 Bulls-Eye for Athletes: Rasch analysis

These packages are used for all upcoming analyses (for the other instruments as well).

Code
#preparing data frame with the bulls eye for athletes items
BE_6items_data_rasch <- bulls_eye_merged_superfinal_data[, c(
  "BE_Competition",
  "BE_Training",
  "BE_PreparationRecovery",
  "BE_Life",
  "BE_Obstacles_Sport",
  "BE_Obstacles_Life"
)]
Code
#investigate missing data
RImissing(BE_6items_data_rasch)

Code
#the following means that participants have answered at least one item
filter <- dplyr::filter

min.responses <- 1

BE_6items_data_rasch <- BE_6items_data_rasch %>%
  filter(6 - rowSums(is.na(
    select(
      .,
      "BE_Competition",
      "BE_Training",
      "BE_PreparationRecovery",
      "BE_Life",
      "BE_Obstacles_Sport",
      "BE_Obstacles_Life"
    )
  )) >= min.responses)
Code
#tile plot (response category frequency for the items)
RItileplot(BE_6items_data_rasch)

A test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta.

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(BE_6items_data_rasch, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 44.7 35 0.13      
Code
#conditional item fit (assessing unidimensionaliity)
simfit1_BE_6items <- RIgetfit(BE_6items_data_rasch, iterations = 100, cpu = 10) 

RIitemfit(BE_6items_data_rasch, simfit1_BE_6items, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
BE_Competition 1.146 [0.815, 1.214] 1.181 [0.791, 1.388] no misfit no misfit -0.18
BE_Training 1.009 [0.834, 1.205] 1.007 [0.821, 1.253] no misfit no misfit -0.20
BE_PreparationRecovery 0.874 [0.808, 1.201] 0.883 [0.798, 1.245] no misfit no misfit -0.19
BE_Life 1.014 [0.861, 1.224] 1.016 [0.865, 1.243] no misfit no misfit -0.37
BE_Obstacles_Sport 0.972 [0.777, 1.229] 0.948 [0.78, 1.232] no misfit no misfit 0.13
BE_Obstacles_Life 1 [0.8, 1.175] 0.994 [0.814, 1.196] no misfit no misfit -0.42
Note:
MSQ values based on conditional calculations (n = 146 complete cases).
Simulation based thresholds from 84 simulated datasets.
Code
#conditional item characteristic curves 
ICCplot(
  as.data.frame(BE_6items_data_rasch),
  itemnumber = 1:4,
  itemdescrip = c("Competition", "Training", "Preparation & Recovery", "Life"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(BE_6items_data_rasch),
  itemnumber = 5:6,
  itemdescrip = c("Obstacles Sport", "Obstacles Life"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
#PCA of residuals (assessing unidimensionality)
RIpcmPCA(na.omit(BE_6items_data_rasch)) 
Eigenvalues Proportion of variance
1.86 32.1%
1.38 22.5%
1.07 16.5%
0.93 16.1%
0.75 12.6%
Code
#item restscores
RIrestscore(BE_6items_data_rasch)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
BE_Competition 0.33 0.33 0.00 0.992 0.02 -0.18
BE_Training 0.29 0.29 0.00 0.992 0.01 -0.20
BE_PreparationRecovery 0.38 0.30 0.08 0.992 0.02 -0.19
BE_Life 0.30 0.31 0.01 0.992 -0.16 -0.37
BE_Obstacles_Sport 0.33 0.31 0.02 0.992 0.33 0.13
BE_Obstacles_Life 0.32 0.32 0.00 0.992 -0.22 -0.42
Code
#loadings on first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(na.omit(BE_6items_data_rasch))

Code
#residual correlations
resid_cor_BE_6items <- RIgetResidCor(BE_6items_data_rasch, iterations = 800, cpu = 10)

RIresidcorr(BE_6items_data_rasch, cutoff = resid_cor_BE_6items$p99) 
BE_Competition BE_Training BE_PreparationRecovery BE_Life BE_Obstacles_Sport BE_Obstacles_Life
BE_Competition
BE_Training 0.09
BE_PreparationRecovery -0.08 -0.11
BE_Life -0.22 -0.11 -0.16
BE_Obstacles_Sport -0.23 -0.25 0 -0.22
BE_Obstacles_Life -0.44 -0.26 -0.09 0.04 0.12
Note:
Relative cut-off value is 0.138, which is 0.267 above the average correlation (-0.129).
Correlations above the cut-off are highlighted in red text.
Code
#partial gamma coefficients (assessing local independence)

clean_names <- janitor::clean_names
filter <- dplyr::filter

RIpartgamLD(BE_6items_data_rasch) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
Code
partgam_LD(as.data.frame(BE_6items_data_rasch)) #shows all output, also non-significant results
                    Item1                  Item2   gamma     se pvalue padj.BH
1          BE_Competition            BE_Training  0.2089 0.1315 0.1121  1.0000
2          BE_Competition BE_PreparationRecovery  0.1289 0.1175 0.2727  1.0000
3          BE_Competition                BE_Life  0.0097 0.1262 0.9387  1.0000
4          BE_Competition     BE_Obstacles_Sport -0.0534 0.1332 0.6882  1.0000
5          BE_Competition      BE_Obstacles_Life -0.2618 0.1222 0.0322  0.9663
6             BE_Training BE_PreparationRecovery -0.0806 0.1123 0.4727  1.0000
7             BE_Training                BE_Life -0.0178 0.1293 0.8905  1.0000
8             BE_Training     BE_Obstacles_Sport -0.1873 0.1265 0.1387  1.0000
9             BE_Training      BE_Obstacles_Life -0.1478 0.1176 0.2089  1.0000
10 BE_PreparationRecovery                BE_Life -0.0254 0.1281 0.8429  1.0000
11 BE_PreparationRecovery     BE_Obstacles_Sport  0.0025 0.1223 0.9839  1.0000
12 BE_PreparationRecovery      BE_Obstacles_Life -0.0724 0.1225 0.5546  1.0000
13                BE_Life     BE_Obstacles_Sport -0.1949 0.1362 0.1525  1.0000
14                BE_Life      BE_Obstacles_Life  0.1802 0.1195 0.1316  1.0000
15     BE_Obstacles_Sport      BE_Obstacles_Life  0.2875 0.1141 0.0117  0.3524
   sig   lower   upper
1      -0.0488  0.4666
2      -0.1014  0.3592
3      -0.2377  0.2571
4      -0.3144  0.2075
5      -0.5013 -0.0222
6      -0.3008  0.1395
7      -0.2713  0.2357
8      -0.4353  0.0606
9      -0.3784  0.0828
10     -0.2764  0.2257
11     -0.2373  0.2422
12     -0.3125  0.1677
13     -0.4620  0.0721
14     -0.0540  0.4145
15      0.0639  0.5112

                    Item1                  Item2   gamma     se pvalue padj.BH
1             BE_Training         BE_Competition  0.1952 0.1150 0.0897  1.0000
2  BE_PreparationRecovery         BE_Competition  0.1799 0.1093 0.0998  1.0000
3  BE_PreparationRecovery            BE_Training  0.0317 0.1119 0.7769  1.0000
4                 BE_Life         BE_Competition  0.0769 0.1051 0.4642  1.0000
5                 BE_Life            BE_Training  0.0201 0.1171 0.8639  1.0000
6                 BE_Life BE_PreparationRecovery -0.1082 0.1190 0.3629  1.0000
7      BE_Obstacles_Sport         BE_Competition -0.1304 0.1179 0.2684  1.0000
8      BE_Obstacles_Sport            BE_Training -0.1291 0.1367 0.3450  1.0000
9      BE_Obstacles_Sport BE_PreparationRecovery  0.1530 0.1321 0.2467  1.0000
10     BE_Obstacles_Sport                BE_Life -0.1940 0.1234 0.1161  1.0000
11      BE_Obstacles_Life         BE_Competition -0.3671 0.1071 0.0006  0.0182
12      BE_Obstacles_Life            BE_Training -0.1755 0.1137 0.1228  1.0000
13      BE_Obstacles_Life BE_PreparationRecovery -0.0439 0.1248 0.7248  1.0000
14      BE_Obstacles_Life                BE_Life  0.1898 0.1130 0.0930  1.0000
15      BE_Obstacles_Life     BE_Obstacles_Sport  0.3054 0.1047 0.0035  0.1061
   sig   lower   upper
1      -0.0302  0.4206
2      -0.0343  0.3941
3      -0.1876  0.2510
4      -0.1291  0.2829
5      -0.2094  0.2495
6      -0.3414  0.1250
7      -0.3614  0.1006
8      -0.3971  0.1389
9      -0.1059  0.4119
10     -0.4359  0.0480
11   * -0.5769 -0.1573
12     -0.3984  0.0474
13     -0.2885  0.2006
14     -0.0316  0.4113
15      0.1002  0.5107
Code
#targeting
RItargeting(BE_6items_data_rasch)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(BE_6items_data_rasch), item = NA)

RIitemHierarchy(BE_6items_data_rasch)

Code
#threshold locations 
RIitemparams(BE_6items_data_rasch)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Threshold 6 Item location
BE_Competition 0.15 -0.58 -0.92 0.45 0.72 0.30 0.02
BE_Training -0.58 -1.91 -0.71 0.55 0.54 2.13 0.01
BE_PreparationRecovery -1.56 -0.84 -0.50 0.21 0.62 2.17 0.02
BE_Life -1.37 -1.24 -0.17 0.09 0.83 0.90 -0.16
BE_Obstacles_Sport -1.02 -1.16 0.40 0.38 1.11 2.29 0.33
BE_Obstacles_Life -2.03 -1.18 0.16 0.42 0.61 0.73 -0.22
Note:
Item location is the average of the thresholds for each item.
Code
#analysis of response categories
RIitemCats(BE_6items_data_rasch, items = "all", xlims = c(-5, 5), legend = "left") 

The first output shows without sample PSI and second output with sample PSI.

Code
#investigating reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(BE_6items_data_rasch))

Code
RItif(na.omit(BE_6items_data_rasch), samplePSI = TRUE)

Based on the analysis of response categories these will be recoded and then undergo Rasch analysis again. Response categories will be recoded with the aim of them being in order, and scale steps will also be merged so that each kept scale step has the highest probability of being used at some point on the logit scale.

2.1 Re-analysis of Bulls Eye recoded version

Based on the analysis of the response categories and threshold locations, scale steps were merged and recoded so that each kept scale step had the highest response probability at some point on the person location scale, and also so that all threshold locations were in the right order for all items. Recoded items were given new names with the abbreviation “_recoded”.

Code
#adding recoded bulls eye items
recode <- dplyr::recode

bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(
    # Recoding BE_Competition (0+1; 2+3; 4+5)
    BE_Competition_recoded = recode(
      BE_Competition,
      `0` = 0, `1` = 0,
      `2` = 1, `3` = 1,
      `4` = 2, `5` = 2,
      `6` = 3
    ),
    # Recoding BE_Training (0+1; 3+4)
    BE_Training_recoded = recode(
      BE_Training,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3,
      `6` = 4
    ),
    #Recoding BE_PreparationRecovery(0+1)
    BE_PreparationRecovery_recoded = recode(
      BE_PreparationRecovery,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    ),
    #Recoding BE_Life (0+1)
    BE_Life_recoded = recode(
      BE_Life,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    ),
    #Recoding BE_Obstacles_Sport (0+1; 3+4; 5+6)
    BE_Obstacles_Sport_recoded = recode(
      BE_Obstacles_Sport,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3, `6` = 3
    ),
    #Recoding BE_Obstacles_Life (0+1)
    BE_Obstacles_Life_recoded = recode(
      BE_Obstacles_Life,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    )
  )
Code
#preparing item data frame for BE recoded

BE_6items_data_rasch_recoded <- bulls_eye_merged_superfinal_data[, c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

min.responses <- 1

BE_6items_data_rasch_recoded <- BE_6items_data_rasch_recoded %>%
  filter(6 - rowSums(is.na(
    select(
      .,
      "BE_Competition_recoded",
      "BE_Training_recoded",
      "BE_PreparationRecovery_recoded",
      "BE_Life_recoded",
      "BE_Obstacles_Sport_recoded",
      "BE_Obstacles_Life_recoded"
    )
  )) >= min.responses) 
Code
#re-analysis of response categories using BE recoded
RIitemCats(BE_6items_data_rasch_recoded, items = "all", xlims = c(-5, 5), legend = "left")

Code
#tile plot (response category frequency for the recoded items)
RItileplot(BE_6items_data_rasch_recoded)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(BE_6items_data_rasch_recoded, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue  sig 
overall 57.9 24 0.00013  ***
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_BE_6items_recoded <- RIgetfit(BE_6items_data_rasch_recoded, iterations = 100, cpu = 10) 

RIitemfit(BE_6items_data_rasch_recoded, simfit1_BE_6items_recoded, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
BE_Competition_recoded 0.97 [0.835, 1.206] 0.966 [0.831, 1.229] no misfit no misfit -0.05
BE_Training_recoded 1.042 [0.822, 1.238] 1.053 [0.801, 1.229] no misfit no misfit 0.03
BE_PreparationRecovery_recoded 0.959 [0.774, 1.213] 0.965 [0.752, 1.243] no misfit no misfit 0.18
BE_Life_recoded 1.076 [0.826, 1.251] 1.078 [0.814, 1.254] no misfit no misfit -0.07
BE_Obstacles_Sport_recoded 0.936 [0.765, 1.242] 0.967 [0.756, 1.296] no misfit no misfit 0.00
BE_Obstacles_Life_recoded 1.005 [0.831, 1.159] 1.021 [0.83, 1.17] no misfit no misfit -0.02
Note:
MSQ values based on conditional calculations (n = 146 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(BE_6items_data_rasch_recoded),
  itemnumber = 1:4,
  itemdescrip = c(
    "Competition_recoded",
    "Training_recoded",
    "Preparation & Recovery_recoded",
    "Life_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(BE_6items_data_rasch_recoded),
  itemnumber = 5:6,
  itemdescrip = c("Obstacles Sport_recoded", "Obstacles Life_recoded"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
#re-analyzing PCA of residuals (assessing unidimensionality)
RIpcmPCA(na.omit(BE_6items_data_rasch_recoded)) 
Eigenvalues Proportion of variance
1.80 30.2%
1.45 24.1%
1.08 17.4%
0.97 16.5%
0.68 11.6%
Code
#re-analyzing item restscores
RIrestscore(BE_6items_data_rasch_recoded)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
BE_Competition_recoded 0.37 0.31 0.06 0.86 -0.07 -0.05
BE_Training_recoded 0.29 0.31 0.02 0.86 0.01 0.03
BE_PreparationRecovery_recoded 0.36 0.35 0.01 0.86 0.16 0.18
BE_Life_recoded 0.31 0.36 0.05 0.86 -0.09 -0.07
BE_Obstacles_Sport_recoded 0.40 0.30 0.10 0.86 -0.02 0.00
BE_Obstacles_Life_recoded 0.34 0.37 0.03 0.86 -0.04 -0.02
Code
#re-analyzing loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(na.omit(BE_6items_data_rasch_recoded))

Code
#re-analysis of residual correlations
resid_cor_BE_6items_recoded <- RIgetResidCor(BE_6items_data_rasch_recoded,
                                             iterations = 800,
                                             cpu = 10)

RIresidcorr(BE_6items_data_rasch_recoded, cutoff = resid_cor_BE_6items_recoded$p99) 
BE_Competition_recoded BE_Training_recoded BE_PreparationRecovery_recoded BE_Life_recoded BE_Obstacles_Sport_recoded BE_Obstacles_Life_recoded
BE_Competition_recoded
BE_Training_recoded 0.2
BE_PreparationRecovery_recoded 0 -0.08
BE_Life_recoded -0.14 -0.14 -0.27
BE_Obstacles_Sport_recoded -0.08 -0.26 -0.04 -0.23
BE_Obstacles_Life_recoded -0.37 -0.25 -0.23 -0.11 0.11
Note:
Relative cut-off value is 0.143, which is 0.269 above the average correlation (-0.126).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(BE_6items_data_rasch_recoded) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
BE_Competition_recoded BE_Training_recoded 0.444 0.141 0.169 0.720 0.047
BE_Training_recoded BE_Competition_recoded 0.419 0.123 0.177 0.661 0.020
BE_Obstacles_Sport_recoded BE_Obstacles_Life_recoded 0.391 0.111 0.174 0.609 0.012
Code
partgam_LD(as.data.frame(BE_6items_data_rasch_recoded)) #shows all output, also non-significant results
                            Item1                          Item2   gamma     se
1          BE_Competition_recoded            BE_Training_recoded  0.4444 0.1407
2          BE_Competition_recoded BE_PreparationRecovery_recoded  0.1895 0.1266
3          BE_Competition_recoded                BE_Life_recoded  0.0435 0.1249
4          BE_Competition_recoded     BE_Obstacles_Sport_recoded -0.0163 0.1471
5          BE_Competition_recoded      BE_Obstacles_Life_recoded -0.2319 0.1338
6             BE_Training_recoded BE_PreparationRecovery_recoded  0.0354 0.1316
7             BE_Training_recoded                BE_Life_recoded  0.0000 0.1391
8             BE_Training_recoded     BE_Obstacles_Sport_recoded -0.1847 0.1606
9             BE_Training_recoded      BE_Obstacles_Life_recoded -0.1848 0.1267
10 BE_PreparationRecovery_recoded                BE_Life_recoded -0.0282 0.1235
11 BE_PreparationRecovery_recoded     BE_Obstacles_Sport_recoded -0.0559 0.1475
12 BE_PreparationRecovery_recoded      BE_Obstacles_Life_recoded -0.1255 0.1218
13                BE_Life_recoded     BE_Obstacles_Sport_recoded -0.2246 0.1427
14                BE_Life_recoded      BE_Obstacles_Life_recoded  0.0979 0.1189
15     BE_Obstacles_Sport_recoded      BE_Obstacles_Life_recoded  0.3915 0.1107
   pvalue padj.BH sig   lower  upper
1  0.0016  0.0474   *  0.1688 0.7201
2  0.1344  1.0000     -0.0586 0.4377
3  0.7278  1.0000     -0.2014 0.2883
4  0.9119  1.0000     -0.3047 0.2721
5  0.0831  1.0000     -0.4941 0.0304
6  0.7882  1.0000     -0.2226 0.2933
7  1.0000  1.0000     -0.2727 0.2727
8  0.2503  1.0000     -0.4995 0.1301
9  0.1449  1.0000     -0.4332 0.0636
10 0.8193  1.0000     -0.2702 0.2138
11 0.7049  1.0000     -0.3450 0.2333
12 0.3027  1.0000     -0.3643 0.1132
13 0.1155  1.0000     -0.5042 0.0551
14 0.4103  1.0000     -0.1351 0.3309
15 0.0004  0.0122   *  0.1745 0.6085

                            Item1                          Item2   gamma     se
1             BE_Training_recoded         BE_Competition_recoded  0.4190 0.1233
2  BE_PreparationRecovery_recoded         BE_Competition_recoded  0.1256 0.1294
3  BE_PreparationRecovery_recoded            BE_Training_recoded  0.0671 0.1341
4                 BE_Life_recoded         BE_Competition_recoded  0.0246 0.1194
5                 BE_Life_recoded            BE_Training_recoded -0.0122 0.1206
6                 BE_Life_recoded BE_PreparationRecovery_recoded -0.1278 0.1234
7      BE_Obstacles_Sport_recoded         BE_Competition_recoded -0.0646 0.1472
8      BE_Obstacles_Sport_recoded            BE_Training_recoded -0.1852 0.1605
9      BE_Obstacles_Sport_recoded BE_PreparationRecovery_recoded  0.0823 0.1416
10     BE_Obstacles_Sport_recoded                BE_Life_recoded -0.1557 0.1396
11      BE_Obstacles_Life_recoded         BE_Competition_recoded -0.4621 0.1084
12      BE_Obstacles_Life_recoded            BE_Training_recoded -0.2538 0.1205
13      BE_Obstacles_Life_recoded BE_PreparationRecovery_recoded -0.1522 0.1196
14      BE_Obstacles_Life_recoded                BE_Life_recoded  0.0776 0.1061
15      BE_Obstacles_Life_recoded     BE_Obstacles_Sport_recoded  0.2752 0.1233
   pvalue padj.BH  sig   lower   upper
1  0.0007  0.0204    *  0.1773  0.6607
2  0.3317  1.0000      -0.1280  0.3793
3  0.6168  1.0000      -0.1957  0.3298
4  0.8365  1.0000      -0.2093  0.2586
5  0.9195  1.0000      -0.2486  0.2242
6  0.3003  1.0000      -0.3697  0.1140
7  0.6608  1.0000      -0.3532  0.2240
8  0.2485  1.0000      -0.4997  0.1293
9  0.5611  1.0000      -0.1952  0.3598
10 0.2648  1.0000      -0.4293  0.1180
11 0.0000  0.0006  *** -0.6747 -0.2496
12 0.0351  1.0000      -0.4900 -0.0177
13 0.2032  1.0000      -0.3867  0.0822
14 0.4648  1.0000      -0.1304  0.2855
15 0.0256  0.7677       0.0336  0.5168
Code
#targeting
RItargeting(BE_6items_data_rasch_recoded)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(BE_6items_data_rasch_recoded), item = NA)

RIitemHierarchy(BE_6items_data_rasch_recoded)

Code
#threshold locations 
RIitemparams(BE_6items_data_rasch_recoded)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
BE_Competition_recoded -1.20 0.04 0.96 NA NA -0.07
BE_Training_recoded -1.66 -1.44 1.17 1.98 NA 0.01
BE_PreparationRecovery_recoded -0.96 -0.78 -0.01 0.47 2.08 0.16
BE_Life_recoded -1.31 -0.46 -0.13 0.67 0.81 -0.09
BE_Obstacles_Sport_recoded -1.09 -0.43 1.45 NA NA -0.02
BE_Obstacles_Life_recoded -1.37 -0.11 0.21 0.46 0.63 -0.04
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(BE_6items_data_rasch_recoded))

Code
RItif(na.omit(BE_6items_data_rasch_recoded), samplePSI = TRUE)

The significant residual correlation between items BE_Competition_recoded and BE_Training_recoded, also support by the significant gamma coefficient between the two, indicate problems with local dependence. Therefore another two re-analyses will be conducted. First without BE_Competition_recoded, and then without BE_Training_recoded.

2.2 Re-analysis of Bulls Eye recoded version, without the BE_Competition_recoded item

Code
#preparing item data frame for BE recoded, without the BE_Competition_recoded item

BE_6items_data_rasch_recoded_trimmed1 <- bulls_eye_merged_superfinal_data[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

min.responses <- 1

BE_6items_data_rasch_recoded_trimmed1 <- BE_6items_data_rasch_recoded_trimmed1 %>%
  filter(6 - rowSums(is.na(
    select(
      .,
      "BE_Training_recoded",
      "BE_PreparationRecovery_recoded",
      "BE_Life_recoded",
      "BE_Obstacles_Sport_recoded",
      "BE_Obstacles_Life_recoded"
    )
  )) >= min.responses) 
Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(BE_6items_data_rasch_recoded_trimmed1, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 36.2 21 0.021   * 
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_BE_6items_recoded_trimmed1 <- RIgetfit(BE_6items_data_rasch_recoded_trimmed1, iterations = 100, cpu = 10) 

RIitemfit(BE_6items_data_rasch_recoded_trimmed1, simfit1_BE_6items_recoded_trimmed1, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
BE_Training_recoded 1.118 [0.809, 1.226] 1.129 [0.821, 1.237] no misfit no misfit 0.02
BE_PreparationRecovery_recoded 0.991 [0.759, 1.225] 0.989 [0.703, 1.256] no misfit no misfit 0.18
BE_Life_recoded 1.065 [0.844, 1.233] 1.046 [0.827, 1.234] no misfit no misfit -0.07
BE_Obstacles_Sport_recoded 0.935 [0.798, 1.198] 0.971 [0.808, 1.197] no misfit no misfit -0.01
BE_Obstacles_Life_recoded 0.861 [0.804, 1.17] 0.889 [0.801, 1.177] no misfit no misfit -0.02
Note:
MSQ values based on conditional calculations (n = 146 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(BE_6items_data_rasch_recoded_trimmed1),
  itemnumber = 1:4,
  itemdescrip = c(
    "Training_recoded",
    "Preparation & Recovery_recoded",
    "Life_recoded",
    "Obstacles Sport_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(BE_6items_data_rasch_recoded_trimmed1),
  itemnumber = 5,
  itemdescrip = c("Obstacles Life_recoded"),
  method = "cut"
)

Code
#re-analyzing PCA of residuals (assessing unidimensionality)
RIpcmPCA(na.omit(BE_6items_data_rasch_recoded_trimmed1)) 
Eigenvalues Proportion of variance
1.57 31.7%
1.48 29.2%
1.03 21.3%
0.90 17.6%
0.01 0.3%
Code
#re-analyzing item restscores
RIrestscore(BE_6items_data_rasch_recoded_trimmed1)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
BE_Training_recoded 0.21 0.30 0.09 0.716 0.00 0.02
BE_PreparationRecovery_recoded 0.32 0.33 0.01 0.855 0.15 0.18
BE_Life_recoded 0.30 0.34 0.04 0.754 -0.09 -0.07
BE_Obstacles_Sport_recoded 0.39 0.29 0.10 0.716 -0.03 -0.01
BE_Obstacles_Life_recoded 0.39 0.34 0.05 0.716 -0.04 -0.02
Code
#re-analyzing loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(na.omit(BE_6items_data_rasch_recoded_trimmed1))

Code
#re-analysis of residual correlations
resid_cor_BE_6items_recoded_trimmed1 <- RIgetResidCor(BE_6items_data_rasch_recoded_trimmed1,
                                             iterations = 800,
                                             cpu = 10)

RIresidcorr(BE_6items_data_rasch_recoded_trimmed1, cutoff = resid_cor_BE_6items_recoded_trimmed1$p99) 
BE_Training_recoded BE_PreparationRecovery_recoded BE_Life_recoded BE_Obstacles_Sport_recoded BE_Obstacles_Life_recoded
BE_Training_recoded
BE_PreparationRecovery_recoded -0.04
BE_Life_recoded -0.11 -0.27
BE_Obstacles_Sport_recoded -0.23 -0.03 -0.24
BE_Obstacles_Life_recoded -0.24 -0.26 -0.16 0.09
Note:
Relative cut-off value is 0.13, which is 0.28 above the average correlation (-0.15).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(BE_6items_data_rasch_recoded_trimmed1) #shows only significant results
[1] "No statistically significant local dependency found."
Code
partgam_LD(as.data.frame(BE_6items_data_rasch_recoded_trimmed1)) #shows all output, also non-significant results
                            Item1                          Item2   gamma     se
1             BE_Training_recoded BE_PreparationRecovery_recoded  0.0824 0.1330
2             BE_Training_recoded                BE_Life_recoded  0.0959 0.1386
3             BE_Training_recoded     BE_Obstacles_Sport_recoded -0.1498 0.1508
4             BE_Training_recoded      BE_Obstacles_Life_recoded -0.2272 0.1238
5  BE_PreparationRecovery_recoded                BE_Life_recoded -0.0766 0.1150
6  BE_PreparationRecovery_recoded     BE_Obstacles_Sport_recoded -0.0516 0.1375
7  BE_PreparationRecovery_recoded      BE_Obstacles_Life_recoded -0.1472 0.1147
8                 BE_Life_recoded     BE_Obstacles_Sport_recoded -0.2977 0.1333
9                 BE_Life_recoded      BE_Obstacles_Life_recoded  0.0280 0.1053
10     BE_Obstacles_Sport_recoded      BE_Obstacles_Life_recoded  0.2998 0.1015
   pvalue padj.BH sig   lower   upper
1  0.5358  1.0000     -0.1783  0.3430
2  0.4890  1.0000     -0.1757  0.3675
3  0.3203  1.0000     -0.4454  0.1457
4  0.0666  1.0000     -0.4699  0.0155
5  0.5052  1.0000     -0.3020  0.1488
6  0.7076  1.0000     -0.3212  0.2180
7  0.1995  1.0000     -0.3720  0.0777
8  0.0255  0.5097     -0.5588 -0.0365
9  0.7900  1.0000     -0.1783  0.2344
10 0.0031  0.0627   .  0.1009  0.4987

                            Item1                          Item2   gamma     se
1  BE_PreparationRecovery_recoded            BE_Training_recoded  0.1033 0.1243
2                 BE_Life_recoded            BE_Training_recoded  0.0966 0.1322
3                 BE_Life_recoded BE_PreparationRecovery_recoded -0.1063 0.1112
4      BE_Obstacles_Sport_recoded            BE_Training_recoded -0.0909 0.1562
5      BE_Obstacles_Sport_recoded BE_PreparationRecovery_recoded  0.1349 0.1274
6      BE_Obstacles_Sport_recoded                BE_Life_recoded -0.1868 0.1347
7       BE_Obstacles_Life_recoded            BE_Training_recoded -0.1429 0.1278
8       BE_Obstacles_Life_recoded BE_PreparationRecovery_recoded -0.0970 0.1177
9       BE_Obstacles_Life_recoded                BE_Life_recoded  0.1588 0.1171
10      BE_Obstacles_Life_recoded     BE_Obstacles_Sport_recoded  0.3301 0.1238
   pvalue padj.BH sig   lower  upper
1  0.4062   1.000     -0.1404 0.3469
2  0.4649   1.000     -0.1625 0.3557
3  0.3392   1.000     -0.3243 0.1117
4  0.5606   1.000     -0.3971 0.2153
5  0.2898   1.000     -0.1149 0.3846
6  0.1656   1.000     -0.4508 0.0773
7  0.2637   1.000     -0.3934 0.1076
8  0.4098   1.000     -0.3277 0.1337
9  0.1748   1.000     -0.0706 0.3883
10 0.0076   0.153      0.0875 0.5727
Code
#targeting
RItargeting(BE_6items_data_rasch_recoded_trimmed1)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(BE_6items_data_rasch_recoded_trimmed1),
                         item = NA)

RIitemHierarchy(BE_6items_data_rasch_recoded_trimmed1)

Code
#threshold locations 
RIitemparams(BE_6items_data_rasch_recoded_trimmed1)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
BE_Training_recoded -1.65 -1.45 1.16 1.95 NA 0
BE_PreparationRecovery_recoded -0.95 -0.78 -0.01 0.45 2.05 0.15
BE_Life_recoded -1.30 -0.46 -0.14 0.66 0.78 -0.09
BE_Obstacles_Sport_recoded -1.08 -0.43 1.43 NA NA -0.03
BE_Obstacles_Life_recoded -1.36 -0.11 0.20 0.44 0.60 -0.04
Note:
Item location is the average of the thresholds for each item.

The first output shows without sample PSI and second output with sample PSI.

Code
#re-analysis reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(BE_6items_data_rasch_recoded_trimmed1))

Code
RItif(na.omit(BE_6items_data_rasch_recoded_trimmed1), samplePSI = TRUE)

2.3 Re-analysis of Bulls Eye recoded version, without the BE_Training_recoded item

Code
#preparing item data frame for BE recoded, without the BE_Training_recoded item
#the BE_Competition_recoded item is put back

BE_6items_data_rasch_recoded_trimmed2 <- bulls_eye_merged_superfinal_data[, c(
  "BE_Competition_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

min.responses <- 1

BE_6items_data_rasch_recoded_trimmed2 <- BE_6items_data_rasch_recoded_trimmed2 %>%
  filter(6 - rowSums(is.na(
    select(
      .,
      "BE_Competition_recoded",
      "BE_PreparationRecovery_recoded",
      "BE_Life_recoded",
      "BE_Obstacles_Sport_recoded",
      "BE_Obstacles_Life_recoded"
    )
  )) >= min.responses) 
Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(BE_6items_data_rasch_recoded_trimmed2, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue  sig 
overall 47.2 20 0.00056  ***
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_BE_6items_recoded_trimmed2 <- RIgetfit(BE_6items_data_rasch_recoded_trimmed2, iterations = 100, cpu = 10) 

RIitemfit(BE_6items_data_rasch_recoded_trimmed2, simfit1_BE_6items_recoded_trimmed2, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
BE_Competition_recoded 1.052 [0.808, 1.18] 1.067 [0.778, 1.2] no misfit no misfit -0.05
BE_PreparationRecovery_recoded 0.985 [0.813, 1.178] 0.989 [0.82, 1.175] no misfit no misfit 0.18
BE_Life_recoded 1.084 [0.785, 1.187] 1.11 [0.751, 1.253] no misfit no misfit -0.06
BE_Obstacles_Sport_recoded 0.904 [0.85, 1.184] 0.922 [0.835, 1.191] no misfit no misfit 0.00
BE_Obstacles_Life_recoded 0.98 [0.776, 1.144] 0.953 [0.76, 1.152] no misfit no misfit -0.01
Note:
MSQ values based on conditional calculations (n = 146 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(BE_6items_data_rasch_recoded_trimmed2),
  itemnumber = 1:4,
  itemdescrip = c(
    "Competition_recoded",
    "Preparation & Recovery_recoded",
    "Life_recoded",
    "Obstacles Sport_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(BE_6items_data_rasch_recoded_trimmed2),
  itemnumber = 5,
  itemdescrip = c("Obstacles Life_recoded"),
  method = "cut"
)

Code
#re-analyzing PCA of residuals (assessing unidimensionality)
RIpcmPCA(na.omit(BE_6items_data_rasch_recoded_trimmed2)) 
Eigenvalues Proportion of variance
1.64 32.7%
1.45 29.1%
1.06 21.4%
0.85 16.6%
0.01 0.1%
Code
#re-analyzing item restscores
RIrestscore(BE_6items_data_rasch_recoded_trimmed2)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
BE_Competition_recoded 0.31 0.32 0.01 0.941 -0.07 -0.05
BE_PreparationRecovery_recoded 0.34 0.35 0.01 0.941 0.16 0.18
BE_Life_recoded 0.32 0.36 0.04 0.941 -0.08 -0.06
BE_Obstacles_Sport_recoded 0.42 0.31 0.11 0.674 -0.02 0.00
BE_Obstacles_Life_recoded 0.36 0.37 0.01 0.941 -0.03 -0.01
Code
#re-analyzing loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(na.omit(BE_6items_data_rasch_recoded_trimmed2))

Code
#re-analysis of residual correlations
resid_cor_BE_6items_recoded_trimmed2 <- RIgetResidCor(BE_6items_data_rasch_recoded_trimmed2,
                                                      iterations = 800,
                                                      cpu = 10)

RIresidcorr(BE_6items_data_rasch_recoded_trimmed2, cutoff = resid_cor_BE_6items_recoded_trimmed2$p99) 
BE_Competition_recoded BE_PreparationRecovery_recoded BE_Life_recoded BE_Obstacles_Sport_recoded BE_Obstacles_Life_recoded
BE_Competition_recoded
BE_PreparationRecovery_recoded 0.02
BE_Life_recoded -0.12 -0.29
BE_Obstacles_Sport_recoded -0.08 -0.07 -0.28
BE_Obstacles_Life_recoded -0.36 -0.26 -0.15 0.07
Note:
Relative cut-off value is 0.125, which is 0.279 above the average correlation (-0.154).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(BE_6items_data_rasch_recoded_trimmed2) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
Code
partgam_LD(as.data.frame(BE_6items_data_rasch_recoded_trimmed2)) #shows all output, also non-significant results
                            Item1                          Item2   gamma     se
1          BE_Competition_recoded BE_PreparationRecovery_recoded  0.2056 0.1113
2          BE_Competition_recoded                BE_Life_recoded  0.0675 0.1323
3          BE_Competition_recoded     BE_Obstacles_Sport_recoded -0.0473 0.1518
4          BE_Competition_recoded      BE_Obstacles_Life_recoded -0.2854 0.1170
5  BE_PreparationRecovery_recoded                BE_Life_recoded -0.0491 0.1105
6  BE_PreparationRecovery_recoded     BE_Obstacles_Sport_recoded -0.0712 0.1361
7  BE_PreparationRecovery_recoded      BE_Obstacles_Life_recoded -0.1128 0.1110
8                 BE_Life_recoded     BE_Obstacles_Sport_recoded -0.2158 0.1436
9                 BE_Life_recoded      BE_Obstacles_Life_recoded  0.0664 0.1148
10     BE_Obstacles_Sport_recoded      BE_Obstacles_Life_recoded  0.3238 0.1092
   pvalue padj.BH sig   lower   upper
1  0.0646  1.0000     -0.0125  0.4238
2  0.6098  1.0000     -0.1918  0.3268
3  0.7552  1.0000     -0.3448  0.2502
4  0.0147  0.2949     -0.5148 -0.0560
5  0.6566  1.0000     -0.2656  0.1674
6  0.6008  1.0000     -0.3381  0.1956
7  0.3095  1.0000     -0.3303  0.1047
8  0.1330  1.0000     -0.4973  0.0657
9  0.5627  1.0000     -0.1585  0.2914
10 0.0030  0.0607   .  0.1097  0.5379

                            Item1                          Item2   gamma     se
1  BE_PreparationRecovery_recoded         BE_Competition_recoded  0.2009 0.1261
2                 BE_Life_recoded         BE_Competition_recoded  0.0784 0.1271
3                 BE_Life_recoded BE_PreparationRecovery_recoded -0.0863 0.1162
4      BE_Obstacles_Sport_recoded         BE_Competition_recoded  0.0085 0.1395
5      BE_Obstacles_Sport_recoded BE_PreparationRecovery_recoded  0.0360 0.1352
6      BE_Obstacles_Sport_recoded                BE_Life_recoded -0.2172 0.1334
7       BE_Obstacles_Life_recoded         BE_Competition_recoded -0.3753 0.1092
8       BE_Obstacles_Life_recoded BE_PreparationRecovery_recoded -0.1641 0.1174
9       BE_Obstacles_Life_recoded                BE_Life_recoded  0.1070 0.1110
10      BE_Obstacles_Life_recoded     BE_Obstacles_Sport_recoded  0.2473 0.1177
   pvalue padj.BH sig   lower   upper
1  0.1111  1.0000     -0.0463  0.4481
2  0.5373  1.0000     -0.1707  0.3276
3  0.4576  1.0000     -0.3141  0.1414
4  0.9514  1.0000     -0.2649  0.2819
5  0.7898  1.0000     -0.2289  0.3010
6  0.1036  1.0000     -0.4787  0.0443
7  0.0006  0.0118   * -0.5893 -0.1613
8  0.1623  1.0000     -0.3942  0.0660
9  0.3350  1.0000     -0.1105  0.3246
10 0.0357  0.7139      0.0165  0.4781
Code
#targeting
RItargeting(BE_6items_data_rasch_recoded_trimmed2)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(BE_6items_data_rasch_recoded_trimmed2),
                         item = NA)

RIitemHierarchy(BE_6items_data_rasch_recoded_trimmed2)

Code
#threshold locations 
RIitemparams(BE_6items_data_rasch_recoded_trimmed2)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
BE_Competition_recoded -1.20 0.03 0.97 NA NA -0.07
BE_PreparationRecovery_recoded -0.96 -0.79 -0.01 0.48 2.11 0.16
BE_Life_recoded -1.32 -0.47 -0.14 0.68 0.83 -0.08
BE_Obstacles_Sport_recoded -1.09 -0.43 1.47 NA NA -0.02
BE_Obstacles_Life_recoded -1.37 -0.12 0.21 0.47 0.66 -0.03
Note:
Item location is the average of the thresholds for each item.

The first output shows without sample PSI and second output with sample PSI.

Code
#re-analysis reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(BE_6items_data_rasch_recoded_trimmed2))

Code
RItif(na.omit(BE_6items_data_rasch_recoded_trimmed2), samplePSI = TRUE)

The problem with local dependence is solved when either removing the BE_Competition_recoded or BE_Training_recoded items from BEA. Since BE_Training_recoded have more scale steps kept (five scale steps, which means more information with increased reliability as a consequence) in comparison to the four scale steps of BE_Competition_recoded, the latter is removed and BE_Training_recoded kept. CLR tests looked slightly better when removing BE_Competition_recoded as well although it was still significant. As a final check let´s run the CLR tests with both items removed.

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
#with both BE_Competition_recoded and BE_Training_recoded items removed
clr_tests(BE_6items_data_rasch_recoded[,-c(1,2)], model = "PCM") 

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 30.8 17 0.021   * 

The CLR tests did not improve when removing both items in comparison to only removing BE_Competition_recoded. Removing the item BE_Competition_recoded is therefore considered the final version of BEA in this Rasch analysis, that will also be used in the upcoming DIF and validity analyses.

2.4 Differential Item Functioning (DIF) and test-retest investigation

DIF investigates invariance, and more specifically whether item functioning varies between subgroups of participants or not. DIF variables to be investigated are age, sex, competition level, sport type (individual/team sport), and high/low scorers on the instrument. All DIF variables will be transformed into binary categories. Test-retest reliability will also be investigated within the framework of DIF. Below is data preparation and analysis for each DIF variable.

Code
#Age DIF variable data preparation

median_age <- median(bulls_eye_merged_superfinal_data$Demografi.Alder)
print(median_age)
[1] 23
Code
#median age is 23 

specific_value <- 23

count_value_median_age <- sum(bulls_eye_merged_superfinal_data$Demografi.Alder == specific_value, na.rm = TRUE)

print(count_value_median_age)
[1] 14
Code
#14 participants is 23 years of age 

range(bulls_eye_merged_superfinal_data$Demografi.Alder)
[1] 15 50
Code
#age range is 15 to 50 years

#age variable is split by median using an inclusive split in which
#participants aged 23 is included in the younger group
#creating age_group variable

bulls_eye_merged_superfinal_data$age_group <- ifelse(bulls_eye_merged_superfinal_data$Demografi.Alder <= median_age,
                                                     "Younger",
                                                     "Older")

dif.age <- factor(bulls_eye_merged_superfinal_data$age_group)

#Age DIF analysis
RIpartgamDIF(BE_6items_data_rasch_recoded_trimmed1, dif.age) 
[1] "No statistically significant DIF found."
Code
#Sex DIF variable data preparation
table(bulls_eye_merged_superfinal_data$Demografi.Kon)

 1  2 
66 90 
Code
#66 females and 90 males 

dif.sex <- factor(bulls_eye_merged_superfinal_data$Demografi.Kon)

#Sex DIF analysis
RIpartgamDIF(BE_6items_data_rasch_recoded_trimmed1, dif.sex)
[1] "No statistically significant DIF found."
Code
#creating binary variable for competition level

#available competition levels in the current sample are
#junior elite, sectional, national, international
#these are categorized as "Lower level competition" (junior elite + sectional)
#and "Higher competition level" (national + international) 

bulls_eye_merged_superfinal_data$comp_level_binary <- ifelse(
  bulls_eye_merged_superfinal_data$comp_level %in% c("junior elite", "sectional"), 
  "Lower competition level", 
  "Higher competition level"
)

table(bulls_eye_merged_superfinal_data$comp_level_binary)

Higher competition level  Lower competition level 
                      93                       63 
Code
#93 higher comp level, 63 lower comp level

dif.comp_level <- factor(bulls_eye_merged_superfinal_data$comp_level_binary)

#Competition level DIF analysis
RIpartgamDIF(BE_6items_data_rasch_recoded_trimmed1, dif.comp_level) 
[1] "No statistically significant DIF found."
Code
#Sport type DIF variable data preparation (individual vs team sport)

table(bulls_eye_merged_superfinal_data$Demografi.Idrotttyp)

  1   2 
 19 137 
Code
#19 individual and 137 team sport athletes

dif.sport_type <- factor(bulls_eye_merged_superfinal_data$Demografi.Idrotttyp)

#Sport type DIF analysis
RIpartgamDIF(BE_6items_data_rasch_recoded_trimmed1, dif.sport_type) 
[1] "No statistically significant DIF found."

In the output, Group 1 is low scorers on BEA, and Group 2 is high scorers on BEA. Groups are split by median.

Code
#DIF analysis for high vs low scorers on the bulls eye instrument (group 1 = low scorers; group 2 = high scorers)
iarm::item_obsexp(PCM(BE_6items_data_rasch_recoded_trimmed1))
Score group 1: 
                               mean obs mean exp std.res sig
BE_Training_recoded            1.761    1.694    0.794      
BE_PreparationRecovery_recoded 1.776    1.684    0.752      
BE_Life_recoded                1.687    1.620    0.557      
BE_Obstacles_Sport_recoded     1.179    1.066    1.263      
BE_Obstacles_Life_recoded      1.284    1.246    0.358      

Score group 2: 
                               mean obs mean exp std.res sig
BE_Training_recoded             2.215    2.272   -0.734     
BE_PreparationRecovery_recoded  3.051    3.114   -0.580     
BE_Life_recoded                 3.215    3.235   -0.172     
BE_Obstacles_Sport_recoded      1.911    1.894    0.208     
BE_Obstacles_Life_recoded       3.278    3.155    0.983     
Code
#Test-retest reliability using DIF, data preparation and analysis 

#bulls eye "test2" equals time point 2 testing in the code below

rename <- dplyr::rename

#creating data frame with bulls eye test 2 items selected
BE_6items_data_rasch_test2 <- bulls_eye_merged_superfinal_data[, c(
  "Användarnamn",
  "BE_Competition_test2",
  "BE_Training_test2",
  "BE_PreparationRecovery_test2",
  "BE_Life_test2",
  "BE_Obstacles_Sport_test2",
  "BE_Obstacles_Life_test2"
)]

columns_to_check <- c(
  "BE_Competition_test2",
  "BE_Training_test2",
  "BE_PreparationRecovery_test2",
  "BE_Life_test2",
  "BE_Obstacles_Sport_test2",
  "BE_Obstacles_Life_test2"
)

#filter rows with at least one non-missing value in the specified columns
BE_6items_data_rasch_test2_filtered <- BE_6items_data_rasch_test2[rowSums(!is.na(BE_6items_data_rasch_test2[columns_to_check])) > 0, ]
#69 participants with a bulls eye time point 2 test

#rename time point 2 items as time point 1 + create new IDs for time point 2 tests 
BE_6items_data_rasch_test2_filtered_renamed_items <- BE_6items_data_rasch_test2_filtered %>%
  mutate(`Användarnamn` = paste0(`Användarnamn`, "_duplicate")) %>%
  rename(
    "BE_Competition" = "BE_Competition_test2",
    "BE_Training" = "BE_Training_test2",
    "BE_PreparationRecovery" = "BE_PreparationRecovery_test2",
    "BE_Life" = "BE_Life_test2",
    "BE_Obstacles_Sport" = "BE_Obstacles_Sport_test2",
    "BE_Obstacles_Life" = "BE_Obstacles_Life_test2"
  )

#adding DIF variable with value 1 for bulls eye time point 2 (time point 2 tests)
BE_6items_data_rasch_test2_filtered_renamed_items$test2_BE <- 1

#now onto preparing data frame for bulls eye test 1
BE_6items_data_rasch_test1 <- bulls_eye_merged_superfinal_data[, c(
  "Användarnamn",
  "BE_Competition",
  "BE_Training",
  "BE_PreparationRecovery",
  "BE_Life",
  "BE_Obstacles_Sport",
  "BE_Obstacles_Life"
)]

#adding DIF variable with value 0 for bulls eye time point 2 (time point 1 tests)
BE_6items_data_rasch_test1$test2_BE <- 0

#merging the two data frames
#in this way time point 2 bulls eye tests can be compared with time point 1
#in which time point 2 tests have been added as "new" individuals in the data (duplicates)
BE_6items_test1_test2_merged <- full_join(BE_6items_data_rasch_test1,
                                          BE_6items_data_rasch_test2_filtered_renamed_items)


#applying the same recoded item structures that came out of the initial BE Rasch analysis (bulls eye recoded version)
BE_6items_test1_test2_merged <- BE_6items_test1_test2_merged %>%
  mutate(
    # Recoding BE_Competition (0+1; 2+3; 4+5)
    BE_Competition_recoded = recode(
      BE_Competition,
      `0` = 0, `1` = 0,
      `2` = 1, `3` = 1,
      `4` = 2, `5` = 2,
      `6` = 3
    ),
    # Recoding BE_Training (0+1; 3+4)
    BE_Training_recoded = recode(
      BE_Training,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3,
      `6` = 4
    ),
    #Recoding BE_PreparationRecovery(0+1)
    BE_PreparationRecovery_recoded = recode(
      BE_PreparationRecovery,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    ),
    #Recoding BE_Life (0+1)
    BE_Life_recoded = recode(
      BE_Life,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    ),
    #Recoding BE_Obstacles_Sport (0+1; 3+4; 5+6)
    BE_Obstacles_Sport_recoded = recode(
      BE_Obstacles_Sport,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3, `6` = 3
    ),
    #Recoding BE_Obstacles_Life (0+1)
    BE_Obstacles_Life_recoded = recode(
      BE_Obstacles_Life,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    )
  )


#test-retest DIF investigation
BE_6items_test1_test2_DIF <- BE_6items_test1_test2_merged[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]


dif.test_retest <- factor(BE_6items_test1_test2_merged$test2_BE)

RIpartgamDIF(BE_6items_test1_test2_DIF, dif.test_retest)
[1] "No statistically significant DIF found."

3 Psychological Flexibility in Sport Scale (PFSS): Rasch analysis and BEA validity investigation

The PFSS instrument is used as part of the validity investigation of BEA. First, PFSS will undergo Rasch analysis. Then a validity analysis of correlating BEA and PFSS thetas will be undertaken. Note that PFSS was not completed by the participants from the intervention study (Hockeystudien), therefore more missing values than in the BEA Rasch analysis are to be expected.

(Note regarding the analysis code: Prior to the publication of PFSS (Johles et al., 2020) the working name of the scale was ADAQ - Accepance, Defusion and Action Questionnaire-sport. Therefore the item names and data frames in the code are called “ADAQ”, but these are the same as the PFSS items)

3.1 Rasch analysis of PFSS

Code
#preparing data frame with PFSS items 
ADAQ_data_rasch <- bulls_eye_merged_superfinal_data[, c("ADAQ.1",
                                                        "ADAQ.2",
                                                        "ADAQ.3",
                                                        "ADAQ.4",
                                                        "ADAQ.5",
                                                        "ADAQ.6",
                                                        "ADAQ.7")]

#investigate missing values
RImissing(ADAQ_data_rasch)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete ADAQ 
ADAQ_data_rasch_na_omit <- ADAQ_data_rasch %>%
  select(ADAQ.1, ADAQ.2, ADAQ.3, ADAQ.4, ADAQ.5, ADAQ.6, ADAQ.7) %>% 
  na.omit()

RImissing(ADAQ_data_rasch_na_omit)
[1] "No missing data."
Code
#tile plot (response category frequency for the items)
RItileplot(ADAQ_data_rasch_na_omit)

Item ADAQ.2 has no responses in scale step zero, and will therefore be merged with scale step 1.

Item ADAQ.2 has no responses in scale step zero, and will therefore be merged with scale step 1.

Code
#merging scale step zero and 1 in ADAQ.2 (creating ADAQ.2_revised) 
bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(ADAQ.2_revised = recode(
    ADAQ.2,
    `0` = 0, `1` = 0,
    `2` = 1,
    `3` = 2,
    `4` = 3,
    `5` = 4,
    `6` = 5
  ))

#preparing data frame with ADAQ items including ADAQ.2_revised
ADAQ_data_rasch <- bulls_eye_merged_superfinal_data[, c("ADAQ.1",
                                                        "ADAQ.2_revised",
                                                        "ADAQ.3",
                                                        "ADAQ.4",
                                                        "ADAQ.5",
                                                        "ADAQ.6",
                                                        "ADAQ.7")]

#removing missing values 
ADAQ_data_rasch_na_omit <- ADAQ_data_rasch %>%
  select(ADAQ.1, ADAQ.2_revised, ADAQ.3, ADAQ.4, ADAQ.5, ADAQ.6, ADAQ.7) %>% 
  na.omit()

RImissing(ADAQ_data_rasch_na_omit)
[1] "No missing data."
Code
RItileplot(ADAQ_data_rasch_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(ADAQ_data_rasch_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 28.8 40 0.91      
Code
#conditional item fit (assessing unidimensionaliity)

simfit1_ADAQ <- RIgetfit(ADAQ_data_rasch_na_omit, iterations = 100, cpu = 10) 

RIitemfit(ADAQ_data_rasch_na_omit, simfit1_ADAQ, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
ADAQ.1 1.254 [0.76, 1.197] 1.296 [0.734, 1.197] 0.057 0.099 -1.22
ADAQ.2_revised 0.737 [0.814, 1.215] 0.751 [0.783, 1.313] 0.077 0.032 -1.17
ADAQ.3 0.9 [0.814, 1.202] 0.916 [0.794, 1.436] no misfit no misfit -0.95
ADAQ.4 1.146 [0.795, 1.332] 1.231 [0.8, 1.377] no misfit no misfit -1.27
ADAQ.5 1.199 [0.711, 1.277] 1.266 [0.751, 1.281] no misfit no misfit -1.14
ADAQ.6 0.745 [0.747, 1.309] 0.73 [0.724, 1.387] 0.002 no misfit -0.87
ADAQ.7 0.975 [0.753, 1.226] 0.972 [0.735, 1.332] no misfit no misfit -1.26
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 54 simulated datasets.
Code
#conditional item characteristic curves
ICCplot(
  as.data.frame(ADAQ_data_rasch_na_omit),
  itemnumber = 1:4,
  itemdescrip = c("ADAQ.1", "ADAQ.2", "ADAQ.3", "ADAQ.4"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(ADAQ_data_rasch_na_omit),
  itemnumber = 5:7,
  itemdescrip = c("ADAQ.5", "ADAQ.6", "ADAQ.7"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
#PCA of residuals (assessing unidimensionality)
RIpcmPCA(ADAQ_data_rasch_na_omit) 
Eigenvalues Proportion of variance
1.71 24.9%
1.39 21.3%
1.21 20.3%
1.01 12.9%
0.87 11.3%
Code
#item restscores
RIrestscore(ADAQ_data_rasch_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
ADAQ.1 0.56 0.61 0.05 0.414 -0.10 -1.22
ADAQ.2_revised 0.73 0.61 0.12 0.048 * -0.05 -1.17
ADAQ.3 0.65 0.62 0.03 0.701 0.18 -0.95
ADAQ.4 0.54 0.62 0.08 0.414 -0.15 -1.27
ADAQ.5 0.58 0.63 0.05 0.509 -0.01 -1.14
ADAQ.6 0.74 0.64 0.10 0.079 . 0.26 -0.87
ADAQ.7 0.65 0.63 0.02 0.701 -0.14 -1.26
Code
#loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(ADAQ_data_rasch_na_omit)

Code
#residual correlations
resid_cor_ADAQ <- RIgetResidCor(ADAQ_data_rasch_na_omit, iterations = 800, cpu = 10)

RIresidcorr(ADAQ_data_rasch_na_omit, cutoff = resid_cor_ADAQ$p99) 
ADAQ.1 ADAQ.2_revised ADAQ.3 ADAQ.4 ADAQ.5 ADAQ.6 ADAQ.7
ADAQ.1
ADAQ.2_revised -0.05
ADAQ.3 -0.26 -0.14
ADAQ.4 -0.23 -0.26 0.1
ADAQ.5 -0.23 0.01 -0.2 -0.08
ADAQ.6 -0.14 -0.08 -0.12 -0.31 -0.25
ADAQ.7 -0.16 -0.11 -0.17 -0.25 -0.36 0.16
Note:
Relative cut-off value is 0.161, which is 0.31 above the average correlation (-0.149).
Correlations above the cut-off are highlighted in red text.
Code
#partial gamma coefficients (assessing local independence)
RIpartgamLD(ADAQ_data_rasch_na_omit) #show only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
ADAQ.6 ADAQ.7 0.504 0.136 0.238 0.77 0.008
Code
partgam_LD(as.data.frame(ADAQ_data_rasch_na_omit)) #show all output, also non-significant results
            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1          ADAQ.1 ADAQ.2_revised  0.0282 0.1726 0.8703  1.0000     -0.3101
2          ADAQ.1         ADAQ.3 -0.0492 0.1709 0.7735  1.0000     -0.3840
3          ADAQ.1         ADAQ.4 -0.0417 0.1591 0.7934  1.0000     -0.3534
4          ADAQ.1         ADAQ.5 -0.0191 0.1579 0.9037  1.0000     -0.3286
5          ADAQ.1         ADAQ.6  0.0222 0.1783 0.9008  1.0000     -0.3273
6          ADAQ.1         ADAQ.7  0.1077 0.1561 0.4903  1.0000     -0.1983
7  ADAQ.2_revised         ADAQ.3 -0.0377 0.2004 0.8506  1.0000     -0.4304
8  ADAQ.2_revised         ADAQ.4  0.0606 0.1836 0.7414  1.0000     -0.2993
9  ADAQ.2_revised         ADAQ.5  0.3566 0.1599 0.0257  1.0000      0.0433
10 ADAQ.2_revised         ADAQ.6  0.1818 0.1649 0.2703  1.0000     -0.1415
11 ADAQ.2_revised         ADAQ.7  0.1562 0.1735 0.3678  1.0000     -0.1838
12         ADAQ.3         ADAQ.4  0.2436 0.1711 0.1546  1.0000     -0.0918
13         ADAQ.3         ADAQ.5 -0.1408 0.1503 0.3488  1.0000     -0.4355
14         ADAQ.3         ADAQ.6 -0.0145 0.1762 0.9345  1.0000     -0.3599
15         ADAQ.3         ADAQ.7 -0.1220 0.1794 0.4967  1.0000     -0.4736
16         ADAQ.4         ADAQ.5  0.0296 0.1511 0.8447  1.0000     -0.2665
17         ADAQ.4         ADAQ.6 -0.3611 0.1513 0.0170  0.7135     -0.6576
18         ADAQ.4         ADAQ.7 -0.2687 0.1607 0.0945  1.0000     -0.5835
19         ADAQ.5         ADAQ.6 -0.1294 0.1451 0.3724  1.0000     -0.4138
20         ADAQ.5         ADAQ.7 -0.3333 0.1656 0.0441  1.0000     -0.6579
21         ADAQ.6         ADAQ.7  0.5041 0.1356 0.0002  0.0084  **  0.2384
     upper
1   0.3664
2   0.2857
3   0.2701
4   0.2904
5   0.3717
6   0.4137
7   0.3549
8   0.4205
9   0.6699
10  0.5051
11  0.4963
12  0.5790
13  0.1538
14  0.3309
15  0.2297
16  0.3257
17 -0.0646
18  0.0462
19  0.1549
20 -0.0088
21  0.7699

            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  ADAQ.2_revised         ADAQ.1  0.1241 0.1719 0.4704  1.0000     -0.2129
2          ADAQ.3         ADAQ.1 -0.0850 0.1727 0.6226  1.0000     -0.4234
3          ADAQ.3 ADAQ.2_revised -0.1163 0.1999 0.5608  1.0000     -0.5081
4          ADAQ.4         ADAQ.1 -0.0909 0.1611 0.5726  1.0000     -0.4067
5          ADAQ.4 ADAQ.2_revised -0.1586 0.1792 0.3762  1.0000     -0.5099
6          ADAQ.4         ADAQ.3  0.3691 0.1550 0.0172  0.7242      0.0653
7          ADAQ.5         ADAQ.1  0.0059 0.1527 0.9691  1.0000     -0.2934
8          ADAQ.5 ADAQ.2_revised  0.2500 0.1550 0.1067  1.0000     -0.0537
9          ADAQ.5         ADAQ.3 -0.2000 0.1542 0.1946  1.0000     -0.5022
10         ADAQ.5         ADAQ.4  0.1226 0.1559 0.4318  1.0000     -0.1831
11         ADAQ.6         ADAQ.1  0.1343 0.1662 0.4188  1.0000     -0.1913
12         ADAQ.6 ADAQ.2_revised  0.1966 0.1785 0.2708  1.0000     -0.1533
13         ADAQ.6         ADAQ.3 -0.0484 0.1940 0.8031  1.0000     -0.4287
14         ADAQ.6         ADAQ.4 -0.2667 0.1764 0.1307  1.0000     -0.6125
15         ADAQ.6         ADAQ.5  0.1765 0.1728 0.3072  1.0000     -0.1622
16         ADAQ.7         ADAQ.1 -0.0440 0.1662 0.7911  1.0000     -0.3697
17         ADAQ.7 ADAQ.2_revised -0.0769 0.1810 0.6708  1.0000     -0.4317
18         ADAQ.7         ADAQ.3 -0.1450 0.1775 0.4139  1.0000     -0.4930
19         ADAQ.7         ADAQ.4 -0.1656 0.1712 0.3335  1.0000     -0.5011
20         ADAQ.7         ADAQ.5 -0.2448 0.1613 0.1291  1.0000     -0.5608
21         ADAQ.7         ADAQ.6  0.4437 0.1470 0.0025  0.1070      0.1555
    upper
1  0.4610
2  0.2534
3  0.2755
4  0.2248
5  0.1927
6  0.6729
7  0.3053
8  0.5537
9  0.1022
10 0.4282
11 0.4600
12 0.5464
13 0.3319
14 0.0792
15 0.5152
16 0.2817
17 0.2778
18 0.2029
19 0.1700
20 0.0713
21 0.7319
Code
#targeting
RItargeting(ADAQ_data_rasch_na_omit)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(ADAQ_data_rasch_na_omit), item = NA)

RIitemHierarchy(ADAQ_data_rasch_na_omit)

Code
#threshold locations 
RIitemparams(ADAQ_data_rasch_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Threshold 6 Item location
ADAQ.1 -3.04 -2.07 -1.30 1.21 1.38 3.25 -0.1
ADAQ.2_revised -2.62 -1.61 0.56 0.45 2.97 NA -0.05
ADAQ.3 -2.08 -1.94 -1.23 1.30 1.07 3.94 0.18
ADAQ.4 -3.26 -1.66 -1.24 0.94 0.90 3.44 -0.15
ADAQ.5 -2.44 -1.69 -0.39 1.20 0.54 2.69 -0.01
ADAQ.6 -1.48 -0.42 -1.18 1.22 0.55 2.87 0.26
ADAQ.7 -2.02 -2.42 -0.01 0.61 0.48 2.53 -0.14
Note:
Item location is the average of the thresholds for each item.
Code
#analysis of response categories
RIitemCats(ADAQ_data_rasch_na_omit, items = "all", xlims = c(-5, 5), legend = "left") 

Code
#analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(ADAQ_data_rasch_na_omit))

Code
RItif(na.omit(ADAQ_data_rasch_na_omit), samplePSI = TRUE)

3.2 Re-analysis of PFSS recoded version

Based on the analysis of the response categories and threshold locations, scale steps were merged and recoded so that each kept scale step had the highest response probability at some point on the person location scale, and also so that all threshold locations were in the right order for all items. Recoded items were given new names with the abbreviation “_recoded”.

Code
#add recoded PFSS items 


bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(
    # Recoding ADAQ.1 (0 + 1)
    ADAQ.1_recoded = recode(
      ADAQ.1,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    ),
    # Recoding ADAQ.2 (0 + 1 + 2; 3 + 4)
    ADAQ.2_recoded = recode(
      ADAQ.2,
      `0` = 0, `1` = 0, `2` = 0,
      `3` = 1, `4` = 1,
      `5` = 2,
      `6` = 3
    ),
    # Recoding ADAQ.3 (0 + 1; 3 + 4)
    ADAQ.3_recoded = recode(
      ADAQ.3,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3,
      `6` = 4
    ),
    
    # Recoding ADAQ.4 (0 + 1)
    ADAQ.4_recoded = recode(
      ADAQ.4,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2,
      `4` = 3,
      `5` = 4,
      `6` = 5
    ),
    
    # Recoding ADAQ.5 (0 + 1; 3 + 4)
    ADAQ.5_recoded = recode(
      ADAQ.5,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3,
      `6` = 4
    ),
    
    # Recoding ADAQ.6 (#ADAQ.6: 0 + 1 + 2; 3 + 4)
    ADAQ.6_recoded = recode(
      ADAQ.6,
      `0` = 0, `1` = 0, `2` = 0,
      `3` = 1, `4` = 1,
      `5` = 2,
      `6` = 3
    ),
    
    # Recoding ADAQ.7 (#ADAQ.7: 0 + 1, 3 + 4)
    ADAQ.7_recoded = recode(
      ADAQ.7,
      `0` = 0, `1` = 0,
      `2` = 1,
      `3` = 2, `4` = 2,
      `5` = 3,
      `6` = 4
    )
  )
Code
#preparing data frame for PFSS recoded 
ADAQ_data_rasch_recoded <- bulls_eye_merged_superfinal_data[, c(
  "ADAQ.1_recoded",
  "ADAQ.2_recoded",
  "ADAQ.3_recoded",
  "ADAQ.4_recoded",
  "ADAQ.5_recoded",
  "ADAQ.6_recoded",
  "ADAQ.7_recoded"
)]

#investigating missing values
RImissing(ADAQ_data_rasch_recoded)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete ADAQ 
ADAQ_data_rasch_recoded_na_omit <- ADAQ_data_rasch_recoded %>%
  select(ADAQ.1_recoded, ADAQ.2_recoded, ADAQ.3_recoded, ADAQ.4_recoded, ADAQ.5_recoded, ADAQ.6_recoded, ADAQ.7_recoded) %>% 
  na.omit()

RImissing(ADAQ_data_rasch_recoded_na_omit)
[1] "No missing data."
Code
#re-analysis of response categories using PFSS recoded
RIitemCats(ADAQ_data_rasch_recoded_na_omit, items = "all", xlims = c(-8, 8), legend = "left")

Code
#tile plot (response category frequency for the recoded items) 
RItileplot(ADAQ_data_rasch_recoded_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(ADAQ_data_rasch_recoded_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig
overall 26  27 0.52      
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_ADAQ_recoded <- RIgetfit(ADAQ_data_rasch_recoded_na_omit, iterations = 100, cpu = 10) 

RIitemfit(ADAQ_data_rasch_recoded_na_omit, simfit1_ADAQ_recoded, cutoff = c(0.005, 0.995)) 
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
ADAQ.1_recoded 1.416 [0.743, 1.239] 1.478 [0.746, 1.299] 0.177 0.179 -0.67
ADAQ.2_recoded 0.716 [0.769, 1.33] 0.708 [0.742, 1.382] 0.053 0.034 -0.37
ADAQ.3_recoded 0.771 [0.713, 1.212] 0.741 [0.716, 1.227] no misfit no misfit -0.57
ADAQ.4_recoded 1.355 [0.699, 1.348] 1.449 [0.726, 1.413] 0.007 0.036 -0.69
ADAQ.5_recoded 1.087 [0.74, 1.275] 1.107 [0.73, 1.342] no misfit no misfit -0.71
ADAQ.6_recoded 0.726 [0.731, 1.236] 0.693 [0.684, 1.683] 0.005 no misfit 0.14
ADAQ.7_recoded 0.962 [0.756, 1.272] 0.958 [0.786, 1.31] no misfit no misfit -0.95
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 96 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(ADAQ_data_rasch_recoded_na_omit),
  itemnumber = 1:4,
  itemdescrip = c(
    "ADAQ.1_recoded",
    "ADAQ.2_recoded",
    "ADAQ.3_recoded",
    "ADAQ.4_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(ADAQ_data_rasch_recoded_na_omit),
  itemnumber = 5:7,
  itemdescrip = c("ADAQ.5_recoded", "ADAQ.6_recoded", "ADAQ.7_recoded"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
#re-analyzing of PCA of residuals (assessing unidimensionality)
RIpcmPCA(ADAQ_data_rasch_recoded_na_omit) 
Eigenvalues Proportion of variance
1.66 26%
1.40 23.5%
1.20 18.7%
1.10 13.3%
0.87 9.7%
Code
#re-analysis of item restscores
RIrestscore(ADAQ_data_rasch_recoded_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
ADAQ.1_recoded 0.54 0.67 0.13 0.032 * -0.08 -0.67
ADAQ.2_recoded 0.78 0.63 0.15 0.003 ** 0.21 -0.37
ADAQ.3_recoded 0.75 0.64 0.11 0.033 * 0.02 -0.57
ADAQ.4_recoded 0.55 0.67 0.12 0.090 . -0.11 -0.69
ADAQ.5_recoded 0.62 0.65 0.03 0.623 -0.12 -0.71
ADAQ.6_recoded 0.76 0.64 0.12 0.032 * 0.72 0.14
ADAQ.7_recoded 0.68 0.65 0.03 0.623 -0.37 -0.95
Code
#re-analysis of loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(ADAQ_data_rasch_recoded_na_omit)

Code
#re-analysis of residual correlations
resid_cor_ADAQ_recoded <- RIgetResidCor(ADAQ_data_rasch_recoded_na_omit,
                                        iterations = 800,
                                        cpu = 10)

RIresidcorr(ADAQ_data_rasch_recoded_na_omit, cutoff = resid_cor_ADAQ_recoded$p99) 
ADAQ.1_recoded ADAQ.2_recoded ADAQ.3_recoded ADAQ.4_recoded ADAQ.5_recoded ADAQ.6_recoded ADAQ.7_recoded
ADAQ.1_recoded
ADAQ.2_recoded -0.16
ADAQ.3_recoded -0.29 0.08
ADAQ.4_recoded -0.28 -0.27 -0.01
ADAQ.5_recoded -0.27 -0.03 -0.11 -0.1
ADAQ.6_recoded -0.2 0.03 0 -0.34 -0.12
ADAQ.7_recoded -0.17 -0.04 -0.14 -0.28 -0.31 0.2
Note:
Relative cut-off value is 0.176, which is 0.309 above the average correlation (-0.134).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(ADAQ_data_rasch_recoded_na_omit) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
ADAQ.6_recoded ADAQ.7_recoded 0.54 0.138 0.269 0.811 0.004
Code
partgam_LD(as.data.frame(ADAQ_data_rasch_recoded_na_omit)) #shows all output, also non-significant results
            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  ADAQ.1_recoded ADAQ.2_recoded -0.0345 0.1860 0.8529  1.0000     -0.3990
2  ADAQ.1_recoded ADAQ.3_recoded -0.2952 0.1828 0.1063  1.0000     -0.6535
3  ADAQ.1_recoded ADAQ.4_recoded -0.0816 0.1491 0.5840  1.0000     -0.3739
4  ADAQ.1_recoded ADAQ.5_recoded  0.0000 0.1461 1.0000  1.0000     -0.2863
5  ADAQ.1_recoded ADAQ.6_recoded -0.0541 0.1993 0.7862  1.0000     -0.4447
6  ADAQ.1_recoded ADAQ.7_recoded  0.1011 0.1483 0.4954  1.0000     -0.1896
7  ADAQ.2_recoded ADAQ.3_recoded  0.3256 0.2527 0.1976  1.0000     -0.1696
8  ADAQ.2_recoded ADAQ.4_recoded  0.0794 0.2173 0.7150  1.0000     -0.3466
9  ADAQ.2_recoded ADAQ.5_recoded  0.3162 0.1934 0.1021  1.0000     -0.0629
10 ADAQ.2_recoded ADAQ.6_recoded  0.2778 0.2247 0.2163  1.0000     -0.1626
11 ADAQ.2_recoded ADAQ.7_recoded  0.2182 0.2095 0.2976  1.0000     -0.1923
12 ADAQ.3_recoded ADAQ.4_recoded  0.1533 0.2054 0.4555  1.0000     -0.2493
13 ADAQ.3_recoded ADAQ.5_recoded -0.0515 0.2015 0.7981  1.0000     -0.4465
14 ADAQ.3_recoded ADAQ.6_recoded -0.0309 0.2559 0.9038  1.0000     -0.5324
15 ADAQ.3_recoded ADAQ.7_recoded -0.0957 0.2121 0.6520  1.0000     -0.5113
16 ADAQ.4_recoded ADAQ.5_recoded  0.0307 0.1596 0.8476  1.0000     -0.2821
17 ADAQ.4_recoded ADAQ.6_recoded -0.3617 0.1684 0.0317  1.0000     -0.6917
18 ADAQ.4_recoded ADAQ.7_recoded -0.2067 0.1697 0.2233  1.0000     -0.5394
19 ADAQ.5_recoded ADAQ.6_recoded -0.1240 0.2048 0.5450  1.0000     -0.5254
20 ADAQ.5_recoded ADAQ.7_recoded -0.3780 0.1647 0.0217  0.9111     -0.7008
21 ADAQ.6_recoded ADAQ.7_recoded  0.5397 0.1383 0.0001  0.0040  **  0.2686
     upper
1   0.3300
2   0.0630
3   0.2106
4   0.2863
5   0.3366
6   0.3918
7   0.8208
8   0.5053
9   0.6954
10  0.7181
11  0.6287
12  0.5558
13  0.3434
14  0.4705
15  0.3200
16  0.3435
17 -0.0317
18  0.1260
19  0.2775
20 -0.0553
21  0.8107

            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  ADAQ.2_recoded ADAQ.1_recoded  0.2319 0.1989 0.2437  1.0000     -0.1580
2  ADAQ.3_recoded ADAQ.1_recoded -0.3097 0.1877 0.0989  1.0000     -0.6776
3  ADAQ.3_recoded ADAQ.2_recoded  0.3253 0.2403 0.1758  1.0000     -0.1456
4  ADAQ.4_recoded ADAQ.1_recoded -0.1337 0.1815 0.4614  1.0000     -0.4894
5  ADAQ.4_recoded ADAQ.2_recoded -0.3200 0.1677 0.0564  1.0000     -0.6487
6  ADAQ.4_recoded ADAQ.3_recoded  0.1875 0.1931 0.3316  1.0000     -0.1910
7  ADAQ.5_recoded ADAQ.1_recoded -0.0508 0.1726 0.7682  1.0000     -0.3890
8  ADAQ.5_recoded ADAQ.2_recoded  0.1579 0.2008 0.4316  1.0000     -0.2356
9  ADAQ.5_recoded ADAQ.3_recoded -0.1376 0.2001 0.4916  1.0000     -0.5297
10 ADAQ.5_recoded ADAQ.4_recoded  0.1088 0.1509 0.4710  1.0000     -0.1870
11 ADAQ.6_recoded ADAQ.1_recoded  0.1971 0.1933 0.3080  1.0000     -0.1818
12 ADAQ.6_recoded ADAQ.2_recoded  0.2039 0.2296 0.3745  1.0000     -0.2461
13 ADAQ.6_recoded ADAQ.3_recoded -0.0787 0.2740 0.7741  1.0000     -0.6158
14 ADAQ.6_recoded ADAQ.4_recoded -0.2222 0.1873 0.2356  1.0000     -0.5894
15 ADAQ.6_recoded ADAQ.5_recoded -0.0392 0.2119 0.8532  1.0000     -0.4545
16 ADAQ.7_recoded ADAQ.1_recoded  0.1638 0.1602 0.3064  1.0000     -0.1501
17 ADAQ.7_recoded ADAQ.2_recoded  0.0492 0.2003 0.8060  1.0000     -0.3433
18 ADAQ.7_recoded ADAQ.3_recoded -0.1667 0.2121 0.4319  1.0000     -0.5823
19 ADAQ.7_recoded ADAQ.4_recoded -0.0296 0.1703 0.8621  1.0000     -0.3634
20 ADAQ.7_recoded ADAQ.5_recoded -0.2893 0.1828 0.1135  1.0000     -0.6475
21 ADAQ.7_recoded ADAQ.6_recoded  0.4468 0.1765 0.0114  0.4773      0.1008
    upper
1  0.6218
2  0.0581
3  0.7962
4  0.2220
5  0.0087
6  0.5660
7  0.2873
8  0.5514
9  0.2545
10 0.4046
11 0.5760
12 0.6539
13 0.4585
14 0.1450
15 0.3760
16 0.4778
17 0.4417
18 0.2490
19 0.3042
20 0.0689
21 0.7928
Code
#re-analysis of targeting
RItargeting(ADAQ_data_rasch_recoded_na_omit)

Code
#re-analysis of item hierarchy
itemlabels <- data.frame(itemnr = names(ADAQ_data_rasch_recoded_na_omit), item = NA)

RIitemHierarchy(ADAQ_data_rasch_recoded_na_omit)

Code
#threshold locations 
RIitemparams(ADAQ_data_rasch_recoded_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
ADAQ.1_recoded -3.39 -2.05 0.82 1.17 3.05 -0.08
ADAQ.2_recoded -2.72 0.64 2.72 NA NA 0.21
ADAQ.3_recoded -2.89 -2.17 1.50 3.65 NA 0.02
ADAQ.4_recoded -2.92 -2.01 0.52 0.65 3.23 -0.11
ADAQ.5_recoded -2.61 -1.28 0.98 2.43 NA -0.12
ADAQ.6_recoded -1.41 0.97 2.61 NA NA 0.72
ADAQ.7_recoded -3.31 -1.10 0.66 2.28 NA -0.37
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(ADAQ_data_rasch_recoded_na_omit))

Code
RItif(na.omit(ADAQ_data_rasch_recoded_na_omit), samplePSI = TRUE)

Removing item ADAQ.6_recoded from PFSS as it does not work well and is potentially redundant. Items ADAQ.1_recoded, ADAQ.2_recoded, and ADAQ.4_recoded does not perform that well neither but is kept for now in the next analysis step. ADAQ.6_recoded is the item that performs worst.

3.3 Re-analysis of PFSS recoded and trimmed version

The ADAQ.6_recoded item is dropped as it did not perform well and is potentially redundant. A re-analysis is conducted with the recoded items but without ADAQ.6_recoded, called the PFFS recoded trimmed version.

Code
#preparing data frame for PFSS recoded and trimmed version

ADAQ_data_rasch_recoded_trimmed <- bulls_eye_merged_superfinal_data[, c(
  "ADAQ.1_recoded",
  "ADAQ.2_recoded",
  "ADAQ.3_recoded",
  "ADAQ.4_recoded",
  "ADAQ.5_recoded",
  "ADAQ.7_recoded"
)]

#investigating missing values
RImissing(ADAQ_data_rasch_recoded_trimmed)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete ADAQ
ADAQ_data_rasch_recoded_trimmed_na_omit <- ADAQ_data_rasch_recoded_trimmed %>%
  select(ADAQ.1_recoded, ADAQ.2_recoded, ADAQ.3_recoded, ADAQ.4_recoded, ADAQ.5_recoded, ADAQ.7_recoded) %>% 
  na.omit()

RImissing(ADAQ_data_rasch_recoded_trimmed_na_omit)
[1] "No missing data."
Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(ADAQ_data_rasch_recoded_trimmed_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig
overall 20  24 0.7       
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_ADAQ_recoded_trimmed <- RIgetfit(ADAQ_data_rasch_recoded_trimmed_na_omit, iterations = 100, cpu = 10) 

RIitemfit(ADAQ_data_rasch_recoded_trimmed_na_omit, simfit1_ADAQ_recoded_trimmed, cutoff = c(0.005, 0.995)) 
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
ADAQ.1_recoded 1.337 [0.737, 1.312] 1.365 [0.74, 1.393] 0.025 no misfit -0.60
ADAQ.2_recoded 0.709 [0.797, 1.186] 0.7 [0.797, 1.242] 0.088 0.097 -0.35
ADAQ.3_recoded 0.746 [0.807, 1.258] 0.72 [0.811, 1.318] 0.061 0.091 -0.49
ADAQ.4_recoded 1.22 [0.741, 1.289] 1.282 [0.709, 1.317] no misfit no misfit -0.62
ADAQ.5_recoded 1.028 [0.755, 1.296] 1.037 [0.737, 1.274] no misfit no misfit -0.64
ADAQ.7_recoded 0.992 [0.745, 1.249] 0.999 [0.749, 1.233] no misfit no misfit -0.88
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 98 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(ADAQ_data_rasch_recoded_trimmed_na_omit),
  itemnumber = 1:4,
  itemdescrip = c(
    "ADAQ.1_recoded",
    "ADAQ.2_recoded",
    "ADAQ.3_recoded",
    "ADAQ.4_recoded"
  )
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(ADAQ_data_rasch_recoded_trimmed_na_omit),
  itemnumber = 5:6,
  itemdescrip = c("ADAQ.5_recoded", "ADAQ.7_recoded")
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
#re-analyzing PCA of residuals (assessing unidimensionality)
RIpcmPCA(ADAQ_data_rasch_recoded_trimmed_na_omit) 
Eigenvalues Proportion of variance
1.49 28.8%
1.41 24.2%
1.21 20.9%
1.05 15.2%
0.83 10.8%
Code
#re-analysis of item restscores
RIrestscore(ADAQ_data_rasch_recoded_trimmed_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
ADAQ.1_recoded 0.53 0.64 0.11 0.066 . 0.01 -0.60
ADAQ.2_recoded 0.77 0.61 0.16 0.001 ** 0.26 -0.35
ADAQ.3_recoded 0.75 0.61 0.14 0.017 * 0.12 -0.49
ADAQ.4_recoded 0.56 0.65 0.09 0.298 -0.01 -0.62
ADAQ.5_recoded 0.62 0.63 0.01 0.901 -0.03 -0.64
ADAQ.7_recoded 0.64 0.62 0.02 0.901 -0.27 -0.88
Code
#re-analysis of loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(ADAQ_data_rasch_recoded_trimmed_na_omit)

Code
#re-analysis of residual correlations
resid_cor_ADAQ_recoded_trimmed <- RIgetResidCor(ADAQ_data_rasch_recoded_trimmed_na_omit,
                                                iterations = 800,
                                                cpu = 10)

RIresidcorr(ADAQ_data_rasch_recoded_trimmed_na_omit, cutoff = resid_cor_ADAQ_recoded_trimmed$p99) 
ADAQ.1_recoded ADAQ.2_recoded ADAQ.3_recoded ADAQ.4_recoded ADAQ.5_recoded ADAQ.7_recoded
ADAQ.1_recoded
ADAQ.2_recoded -0.16
ADAQ.3_recoded -0.3 0.11
ADAQ.4_recoded -0.35 -0.3 -0.02
ADAQ.5_recoded -0.3 -0.02 -0.1 -0.14
ADAQ.7_recoded -0.15 0 -0.09 -0.28 -0.27
Note:
Relative cut-off value is 0.167, which is 0.326 above the average correlation (-0.159).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(ADAQ_data_rasch_recoded_trimmed_na_omit) #shows only significant results
[1] "No statistically significant local dependency found."
Code
partgam_LD(as.data.frame(ADAQ_data_rasch_recoded_trimmed_na_omit)) #shows all output, also non-significant results
            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  ADAQ.1_recoded ADAQ.2_recoded -0.0122 0.1796 0.9459  1.0000     -0.3641
2  ADAQ.1_recoded ADAQ.3_recoded -0.2441 0.1829 0.1820  1.0000     -0.6025
3  ADAQ.1_recoded ADAQ.4_recoded -0.1250 0.1561 0.4232  1.0000     -0.4309
4  ADAQ.1_recoded ADAQ.5_recoded -0.0200 0.1526 0.8957  1.0000     -0.3190
5  ADAQ.1_recoded ADAQ.7_recoded  0.1429 0.1507 0.3432  1.0000     -0.1525
6  ADAQ.2_recoded ADAQ.3_recoded  0.3167 0.1967 0.1075  1.0000     -0.0689
7  ADAQ.2_recoded ADAQ.4_recoded -0.1176 0.1816 0.5171  1.0000     -0.4736
8  ADAQ.2_recoded ADAQ.5_recoded  0.3951 0.1691 0.0195  0.5837      0.0637
9  ADAQ.2_recoded ADAQ.7_recoded  0.3333 0.1725 0.0533  1.0000     -0.0048
10 ADAQ.3_recoded ADAQ.4_recoded  0.3000 0.1714 0.0801  1.0000     -0.0360
11 ADAQ.3_recoded ADAQ.5_recoded -0.0161 0.2069 0.9379  1.0000     -0.4216
12 ADAQ.3_recoded ADAQ.7_recoded  0.0946 0.2005 0.6371  1.0000     -0.2984
13 ADAQ.4_recoded ADAQ.5_recoded -0.1810 0.1731 0.2959  1.0000     -0.5202
14 ADAQ.4_recoded ADAQ.7_recoded -0.2444 0.1684 0.1466  1.0000     -0.5745
15 ADAQ.5_recoded ADAQ.7_recoded -0.2077 0.1775 0.2420  1.0000     -0.5557
    upper
1  0.3398
2  0.1144
3  0.1809
4  0.2790
5  0.4383
6  0.7022
7  0.2383
8  0.7264
9  0.6715
10 0.6360
11 0.3893
12 0.4876
13 0.1583
14 0.0856
15 0.1402

            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  ADAQ.2_recoded ADAQ.1_recoded  0.3168 0.1727 0.0666  1.0000     -0.0217
2  ADAQ.3_recoded ADAQ.1_recoded -0.1719 0.2070 0.4064  1.0000     -0.5777
3  ADAQ.3_recoded ADAQ.2_recoded  0.3274 0.1995 0.1007  1.0000     -0.0635
4  ADAQ.4_recoded ADAQ.1_recoded -0.0977 0.1660 0.5563  1.0000     -0.4231
5  ADAQ.4_recoded ADAQ.2_recoded -0.4356 0.1512 0.0040  0.1186     -0.7318
6  ADAQ.4_recoded ADAQ.3_recoded  0.0350 0.1938 0.8568  1.0000     -0.3449
7  ADAQ.5_recoded ADAQ.1_recoded -0.0802 0.1496 0.5919  1.0000     -0.3733
8  ADAQ.5_recoded ADAQ.2_recoded  0.1975 0.1829 0.2804  1.0000     -0.1611
9  ADAQ.5_recoded ADAQ.3_recoded -0.1603 0.1908 0.4007  1.0000     -0.5342
10 ADAQ.5_recoded ADAQ.4_recoded -0.0239 0.1506 0.8737  1.0000     -0.3190
11 ADAQ.7_recoded ADAQ.1_recoded  0.1152 0.1574 0.4642  1.0000     -0.1933
12 ADAQ.7_recoded ADAQ.2_recoded  0.1818 0.1874 0.3321  1.0000     -0.1856
13 ADAQ.7_recoded ADAQ.3_recoded  0.0483 0.2191 0.8256  1.0000     -0.3812
14 ADAQ.7_recoded ADAQ.4_recoded -0.0101 0.1712 0.9529  1.0000     -0.3456
15 ADAQ.7_recoded ADAQ.5_recoded -0.0909 0.1986 0.6471  1.0000     -0.4801
     upper
1   0.6552
2   0.2339
3   0.7184
4   0.2277
5  -0.1393
6   0.4148
7   0.2129
8   0.5560
9   0.2136
10  0.2712
11  0.4237
12  0.5492
13  0.4778
14  0.3254
15  0.2983
Code
#re-analysis of targeting
RItargeting(ADAQ_data_rasch_recoded_trimmed_na_omit)

Code
#re-analysis of item hierarchy
itemlabels <- data.frame(itemnr = names(ADAQ_data_rasch_recoded_trimmed_na_omit),
                         item = NA)

RIitemHierarchy(ADAQ_data_rasch_recoded_trimmed_na_omit)

Code
#threshold locations 
RIitemparams(ADAQ_data_rasch_recoded_trimmed_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
ADAQ.1_recoded -3.01 -1.86 0.84 1.12 2.95 0.01
ADAQ.2_recoded -2.49 0.65 2.64 NA NA 0.26
ADAQ.3_recoded -2.51 -2.02 1.47 3.53 NA 0.12
ADAQ.4_recoded -2.55 -1.82 0.54 0.62 3.13 -0.01
ADAQ.5_recoded -2.27 -1.17 0.96 2.35 NA -0.03
ADAQ.7_recoded -2.96 -0.99 0.65 2.21 NA -0.27
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(ADAQ_data_rasch_recoded_trimmed_na_omit))

Code
RItif(na.omit(ADAQ_data_rasch_recoded_trimmed_na_omit), samplePSI = TRUE)

Items ADAQ.1_recoded, ADAQ.2_recoded, and ADAQ.3_recoded still does not perform optimal (see conditional item fit and item restscores) but is kept in the PFSS scale for the validity analysis.

3.4 PFSS and BEA validity analysis

Construct validity between PFSS and BEA is investigated by correlation analysis. As both scales has undergone Rasch analysis, participants theta estimates (computed with weighted likelihood estimation - WLE) will be correlated, and since thetas are interval data Pearson’s correlation is applied.

The following data frame will be used for all validity analyses.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item data frames
BE_validity_analyses <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

ADAQ_validity_analyses <- data_complete_cases[, c(
  "ADAQ.1_recoded",
  "ADAQ.2_recoded",
  "ADAQ.3_recoded",
  "ADAQ.4_recoded",
  "ADAQ.5_recoded",
  "ADAQ.7_recoded"
)]

#extracting thetas
BE_thetas <- RIestThetas(BE_validity_analyses)

ADAQ_thetas <- RIestThetas(ADAQ_validity_analyses)

#adding thetas to the data frame
data_complete_cases$BE_thetas_WLE <- BE_thetas$WLE
data_complete_cases$ADAQ_thetas_WLE <- ADAQ_thetas$WLE
Code
correlation_BE_ADAQ <- cor_test(
  data_complete_cases,
  x = "BE_thetas_WLE",
  y = "ADAQ_thetas_WLE",
  method = "pearson"
)

print(correlation_BE_ADAQ)
Parameter1    |      Parameter2 |    r |       95% CI | t(112) |         p
--------------------------------------------------------------------------
BE_thetas_WLE | ADAQ_thetas_WLE | 0.30 | [0.13, 0.46] |   3.38 | < .001***

Observations: 114
Code
plot(correlation_BE_ADAQ)

Some participants did not complete BEA and PFSS during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#time latency variables have previously been created
#using data frame from correlation calculation to account for removed NA
mean(data_complete_cases$time_latency_survey1_bulls_eye1_days,
     na.rm = TRUE)
[1] 0.877088
Code
sd(data_complete_cases$time_latency_survey1_bulls_eye1_days,
   na.rm = TRUE)
[1] 2.824687
Code
range(data_complete_cases$time_latency_survey1_bulls_eye1_days,
      na.rm = TRUE)
[1]  0.003171296 21.032523148
Code
#count participants with a time latency of more than 1 day 
count_greater_than_1 <- sum(data_complete_cases$time_latency_survey1_bulls_eye1_days > 1,
                            na.rm = TRUE)

#count participants with a time latency of more than 3 days 
count_greater_than_3 <- sum(data_complete_cases$time_latency_survey1_bulls_eye1_days > 3,
                            na.rm = TRUE)

print(paste("Number of participants with a time latency of more than 1 day:", count_greater_than_1))
[1] "Number of participants with a time latency of more than 1 day: 15"
Code
print(paste("Number of participants with a time latency of more than 3 days:", count_greater_than_3))
[1] "Number of participants with a time latency of more than 3 days: 10"
Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#sensitivity analysis ADAQ and BE
BE_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

ADAQ_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "ADAQ.1_recoded",
  "ADAQ.2_recoded",
  "ADAQ.3_recoded",
  "ADAQ.4_recoded",
  "ADAQ.5_recoded",
  "ADAQ.7_recoded"
)]

#extracting thetas
BE_thetas_less1day <- RIestThetas(BE_validity_analyses_less1day)

ADAQ_thetas_less1day <- RIestThetas(ADAQ_validity_analyses_less1day)

#adding thetas to the data frame
data_complete_cases_less_than_1day$BE_thetas_less1day_WLE <- BE_thetas_less1day$WLE
data_complete_cases_less_than_1day$ADAQ_thetas_less1day_WLE <- ADAQ_thetas_less1day$WLE

#analysis only with participants completing Bulls Eye and ADAQ the same day
correlation_BE_ADAQ_less1day <- cor_test(
  data_complete_cases_less_than_1day,
  x = "BE_thetas_less1day_WLE",
  y = "ADAQ_thetas_less1day_WLE",
  method = "pearson"
)

print(correlation_BE_ADAQ_less1day)
Parameter1             |               Parameter2 |    r |       95% CI | t(97) |       p
-----------------------------------------------------------------------------------------
BE_thetas_less1day_WLE | ADAQ_thetas_less1day_WLE | 0.32 | [0.13, 0.49] |  3.36 | 0.001**

Observations: 99
Code
#we can also check if there is a statistically significant difference 
#by adding same or different day as a predictor 
data_complete_cases_sensitivity_analysis <- data_complete_cases %>% 
  mutate(group = factor(ifelse(time_latency_survey1_bulls_eye1_days < 1,"same_day","different_day")))

m0_adaq <- lm(ADAQ_thetas_WLE ~ BE_thetas_WLE, data = data_complete_cases_sensitivity_analysis)
m1_adaq <- lm(ADAQ_thetas_WLE ~ BE_thetas_WLE + group, data = data_complete_cases_sensitivity_analysis)

anova(m0_adaq, m1_adaq)
Analysis of Variance Table

Model 1: ADAQ_thetas_WLE ~ BE_thetas_WLE
Model 2: ADAQ_thetas_WLE ~ BE_thetas_WLE + group
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    112 278.88                           
2    111 275.47  1     3.409 1.3736 0.2437

Adding whether BE and PFSS were completed during the same day or not as a predictor did not significantly alter the results.

4 Satisfaction With Life Scale (SWLS): Rasch analysis and BEA validity investigation

The SWLS instrument is used as part of the validity investigation of BEA. First, SWLS will undergo Rasch analysis. Then a validity analysis of correlating BEA and SWLS thetas will be undertaken. Note that SWLS was not completed by the participants from the intervention study (Hockeystudien), therefore more missing values than in the BEA Rasch analysis are to be expected.

4.1 Rasch analysis of SWLS

Code
#preparing data frame with SWLS items 
SWLS_data_rasch <- bulls_eye_merged_superfinal_data[, c("SWLS.1", "SWLS.2", "SWLS.3", "SWLS.4", "SWLS.5")]


#investigate missing values
RImissing(SWLS_data_rasch)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete SWLS
SWLS_data_rasch_na_omit <- SWLS_data_rasch %>%
  select(SWLS.1, SWLS.2, SWLS.3, SWLS.4, SWLS.5) %>% 
  na.omit()

RImissing(SWLS_data_rasch_na_omit)
[1] "No missing data."
Code
#tile plot (response category frequency for the items)
RItileplot(SWLS_data_rasch_na_omit)

SWLS.2 has no responses in scale step one and will therefore be merged with scale step zero. SWLS.3 and SWLS.4 has no responses in scale step zero and will therefore be merged with scale step one.

SWLS.2 has no responses in scale step one and will therefore be merged with scale step zero. SWLS.3 and SWLS.4 has no responses in scale step zero and will therefore be merged with scale step one.

Code
#merging scale step zero and one in SWLS.2, SWLS.3, and SWLS.4 (creating SWLS.2_revised, SWLS.3_revised, and SWLS.4_revised)
bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(
    # Recoding SWLS.2 (0+1)
    SWLS.2_revised = recode(SWLS.2,
                            `0` = 0, `1` = 0, 
                            `2` = 1,
                            `3` = 2, 
                            `4` = 3,
                            `5` = 4,
                            `6` = 5),
    # Recoding SWLS.3 (0+1)
    SWLS.3_revised = recode(SWLS.3,
                            `0` = 0, `1` = 0, 
                            `2` = 1,
                            `3` = 2, 
                            `4` = 3,
                            `5` = 4,
                            `6` = 5),
    # Recoding SWLS.4 (0+1)
    SWLS.4_revised = recode(SWLS.4,
                            `0` = 0, `1` = 0, 
                            `2` = 1,
                            `3` = 2, 
                            `4` = 3,
                            `5` = 4,
                            `6` = 5))


#preparing data frame with SWLS items including SWLS.2_revised, SWLS.3_revised, and SWLS.4_revised
SWLS_data_rasch <- bulls_eye_merged_superfinal_data[, c("SWLS.1", "SWLS.2_revised", "SWLS.3_revised", "SWLS.4_revised", "SWLS.5")]


#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete SWLS
SWLS_data_rasch_na_omit <- SWLS_data_rasch %>%
  select(SWLS.1, SWLS.2_revised, SWLS.3_revised, SWLS.4_revised, SWLS.5) %>% 
  na.omit()

RImissing(SWLS_data_rasch_na_omit)
[1] "No missing data."
Code
RItileplot(SWLS_data_rasch_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(SWLS_data_rasch_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 15.8 26 0.94      
Code
#conditional item fit (assessing unidimensionaliity)
simfit1_SWLS <- RIgetfit(SWLS_data_rasch_na_omit, iterations = 100, cpu = 10) 

RIitemfit(SWLS_data_rasch_na_omit, simfit1_SWLS, cutoff = c(0.005, 0.995)) 
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
SWLS.1 1.096 [0.819, 1.242] 1.213 [0.797, 1.237] no misfit no misfit -0.48
SWLS.2_revised 1.128 [0.799, 1.211] 1.188 [0.804, 1.308] no misfit no misfit -0.87
SWLS.3_revised 0.814 [0.755, 1.205] 0.714 [0.732, 1.319] no misfit 0.018 -0.86
SWLS.4_revised 0.775 [0.8, 1.272] 0.839 [0.789, 1.292] 0.025 no misfit -0.93
SWLS.5 1.2 [0.828, 1.31] 1.162 [0.769, 1.392] no misfit no misfit -0.32
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 74 simulated datasets.
Code
#conditional item characteristic curves
ICCplot(
  as.data.frame(SWLS_data_rasch_na_omit),
  itemnumber = 1:4,
  itemdescrip = c("SWLS.1", "SWLS.2", "SWLS.3", "SWLS.4"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(SWLS_data_rasch_na_omit),
  itemnumber = 5,
  itemdescrip = c("SWLS.5"),
  method = "cut"
)

Code
#PCA of residuals (assessing unidimensionality)
RIpcmPCA(SWLS_data_rasch_na_omit)
Eigenvalues Proportion of variance
1.49 30.4%
1.34 29.5%
1.27 26%
0.89 14%
0.01 0.1%
Code
#item restscores
RIrestscore(SWLS_data_rasch_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
SWLS.1 0.44 0.50 0.06 0.700 0.19 -0.48
SWLS.2_revised 0.47 0.50 0.03 0.749 -0.20 -0.87
SWLS.3_revised 0.64 0.50 0.14 0.083 . -0.19 -0.86
SWLS.4_revised 0.57 0.49 0.08 0.615 -0.26 -0.93
SWLS.5 0.50 0.52 0.02 0.749 0.35 -0.32
Code
#loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(SWLS_data_rasch_na_omit)

Code
#residual correlations
resid_cor_SWLS <- RIgetResidCor(SWLS_data_rasch_na_omit,
                                iterations = 800,
                                cpu = 10)

RIresidcorr(SWLS_data_rasch_na_omit, cutoff = resid_cor_SWLS$p99) 
SWLS.1 SWLS.2_revised SWLS.3_revised SWLS.4_revised SWLS.5
SWLS.1
SWLS.2_revised -0.24
SWLS.3_revised -0.12 -0.23
SWLS.4_revised -0.17 -0.13 0.12
SWLS.5 -0.33 -0.26 -0.28 -0.28
Note:
Relative cut-off value is 0.127, which is 0.319 above the average correlation (-0.192).
Correlations above the cut-off are highlighted in red text.
Code
#partial gamma coefficients (assessing local independence)
RIpartgamLD(SWLS_data_rasch_na_omit) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
SWLS.3_revised SWLS.4_revised 0.47 0.149 0.178 0.762 0.032
Code
partgam_LD(as.data.frame(SWLS_data_rasch_na_omit)) #shows all output, also non-significant results
            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1          SWLS.1 SWLS.2_revised -0.0373 0.1498 0.8032  1.0000     -0.3310
2          SWLS.1 SWLS.3_revised -0.2050 0.1631 0.2088  1.0000     -0.5247
3          SWLS.1 SWLS.4_revised -0.1682 0.1411 0.2333  1.0000     -0.4449
4          SWLS.1         SWLS.5 -0.1139 0.1370 0.4057  1.0000     -0.3825
5  SWLS.2_revised SWLS.3_revised -0.1416 0.1678 0.3988  1.0000     -0.4705
6  SWLS.2_revised SWLS.4_revised -0.0448 0.1521 0.7684  1.0000     -0.3429
7  SWLS.2_revised         SWLS.5  0.0615 0.1384 0.6568  1.0000     -0.2097
8  SWLS.3_revised SWLS.4_revised  0.4703 0.1489 0.0016  0.0318   *  0.1784
9  SWLS.3_revised         SWLS.5  0.1808 0.1491 0.2252  1.0000     -0.1114
10 SWLS.4_revised         SWLS.5 -0.0563 0.1498 0.7068  1.0000     -0.3499
    upper
1  0.2563
2  0.1147
3  0.1084
4  0.1546
5  0.1873
6  0.2533
7  0.3327
8  0.7622
9  0.4730
10 0.2372

            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  SWLS.2_revised         SWLS.1 -0.0120 0.1343 0.9285       1     -0.2753
2  SWLS.3_revised         SWLS.1  0.0588 0.1791 0.7426       1     -0.2922
3  SWLS.3_revised SWLS.2_revised -0.0652 0.1765 0.7118       1     -0.4112
4  SWLS.4_revised         SWLS.1 -0.1897 0.1686 0.2606       1     -0.5201
5  SWLS.4_revised SWLS.2_revised  0.0093 0.1615 0.9543       1     -0.3074
6  SWLS.4_revised SWLS.3_revised  0.2511 0.1606 0.1180       1     -0.0637
7          SWLS.5         SWLS.1  0.0574 0.1271 0.6518       1     -0.1918
8          SWLS.5 SWLS.2_revised  0.0871 0.1479 0.5559       1     -0.2028
9          SWLS.5 SWLS.3_revised  0.1250 0.1575 0.4275       1     -0.1838
10         SWLS.5 SWLS.4_revised -0.0968 0.1484 0.5144       1     -0.3877
    upper
1  0.2512
2  0.4099
3  0.2807
4  0.1407
5  0.3259
6  0.5659
7  0.3066
8  0.3771
9  0.4338
10 0.1941
Code
#targeting
RItargeting(SWLS_data_rasch_na_omit)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(SWLS_data_rasch_na_omit), item = NA)

RIitemHierarchy(SWLS_data_rasch_na_omit) 

Code
#threshold locations 
RIitemparams(SWLS_data_rasch_na_omit) 
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Threshold 6 Item location
SWLS.1 -1.81 -1.13 0.40 -0.59 1.22 3.05 0.19
SWLS.2_revised -2.27 -0.01 -0.90 0.61 1.57 NA -0.2
SWLS.3_revised -1.79 0.47 -1.46 0.33 1.52 NA -0.19
SWLS.4_revised -1.81 -0.71 -0.89 0.54 1.55 NA -0.26
SWLS.5 -0.23 -0.80 0.58 -0.35 0.83 2.06 0.35
Note:
Item location is the average of the thresholds for each item.
Code
#analysis of response categories
RIitemCats(SWLS_data_rasch_na_omit, items = "all", xlims = c(-5, 5), legend = "left") 

Code
#analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(SWLS_data_rasch_na_omit))

Code
RItif(na.omit(SWLS_data_rasch_na_omit), samplePSI = TRUE)

4.2 Re-analysis of SWLS recoded version

Based on the analysis of the response categories and threshold locations, scale steps were merged and recoded so that each kept scale step had the highest response probability at some point on the person location scale, and also so that all threshold locations were in the right order for all items. Recoded items were given new names with the abbrevation “_recoded”.

Code
#add recoded SWLS items 

bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(
    # Recoding SWLS.1 (0 + 1; 2 + 3)
    SWLS.1_recoded = recode(SWLS.1,
                            `0` = 0, `1` = 0,
                            `2` = 1, `3` = 1,
                            `4` = 2,
                            `5` = 3,
                            `6` = 4),
    # Recoding SWLS.2 (0+1+2; 3+4)
    SWLS.2_recoded = recode(SWLS.2,
                            `0` = 0, `1` = 0, `2` = 0,
                            `3` = 1, `4` = 1,
                            `5` = 2,
                            `6` = 3),
    # Recoding SWLS.3 (0+1+2; 3+4)
    SWLS.3_recoded = recode(SWLS.3,
                            `0` = 0, `1` = 0, `2` = 0,
                            `3` = 1, `4` = 1,
                            `5` = 2,
                            `6` = 3),
    # Recoding SWLS.4 (0+1+2+3)
    SWLS.4_recoded = recode(SWLS.4,
                            `0` = 0, `1` = 0, `2` = 0, `3` = 0,
                            `4` = 1,
                            `5` = 2,
                            `6` = 3),
    # Recoding SWLS.5 (0+1, 2+3)
    SWLS.5_recoded = recode(SWLS.5,
                            `0` = 0, `1` = 0,
                            `2` = 1, `3` = 1,
                            `4` = 2,
                            `5` = 3,
                            `6` = 4))
Code
#preparing data frame for SWLS recoded
SWLS_data_rasch_recoded <- bulls_eye_merged_superfinal_data[, c(
  "SWLS.1_recoded",
  "SWLS.2_recoded",
  "SWLS.3_recoded",
  "SWLS.4_recoded",
  "SWLS.5_recoded"
)]

#investigate missing values 
RImissing(SWLS_data_rasch_recoded)

Code
#removing missing values
SWLS_data_rasch_recoded_na_omit <- SWLS_data_rasch_recoded %>%
  select(SWLS.1_recoded, SWLS.2_recoded, SWLS.3_recoded, SWLS.4_recoded, SWLS.5_recoded) %>% 
  na.omit()

RImissing(SWLS_data_rasch_recoded_na_omit)
[1] "No missing data."
Code
#re-analysis of response categories using SWLS recoded
RIitemCats(SWLS_data_rasch_recoded_na_omit, items = "all", xlims = c(-8, 8), legend = "left") 

Code
#tile plot (response category frequency for the recoded items) 
RItileplot(SWLS_data_rasch_recoded_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(SWLS_data_rasch_recoded_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 14.4 16 0.57      
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_SWLS_recoded <- RIgetfit(SWLS_data_rasch_recoded_na_omit, iterations = 100, cpu = 10) 

RIitemfit(SWLS_data_rasch_recoded_na_omit, simfit1_SWLS_recoded, cutoff = c(0.005, 0.995)) 
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
SWLS.1_recoded 1.146 [0.793, 1.277] 1.154 [0.796, 1.265] no misfit no misfit 0.14
SWLS.2_recoded 1.09 [0.793, 1.237] 1.173 [0.8, 1.269] no misfit no misfit -0.29
SWLS.3_recoded 0.769 [0.769, 1.219] 0.725 [0.729, 1.346] no misfit 0.004 -0.40
SWLS.4_recoded 0.851 [0.761, 1.246] 0.853 [0.732, 1.348] no misfit no misfit -0.16
SWLS.5_recoded 1.122 [0.793, 1.194] 1.186 [0.753, 1.203] no misfit no misfit 0.00
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(SWLS_data_rasch_recoded_na_omit),
  itemnumber = 1:4,
  itemdescrip = c(
    "SWLS.1_recoded",
    "SWLS.2_recoded",
    "SWLS.3_recoded",
    "SWLS.4_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(SWLS_data_rasch_recoded_na_omit),
  itemnumber = 5,
  itemdescrip = c("SWLS.5_recoded"),
  method = "cut"
)

Code
#re-analysis of PCA of residuals (assessing unidimensionality)
RIpcmPCA(SWLS_data_rasch_recoded_na_omit)
Eigenvalues Proportion of variance
1.47 30.1%
1.33 29.2%
1.28 25.7%
0.91 14.8%
0.01 0.1%
Code
#re-analysis of item restscores
RIrestscore(SWLS_data_rasch_recoded_na_omit) 
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
SWLS.1_recoded 0.44 0.52 0.08 0.441 0.25 0.14
SWLS.2_recoded 0.48 0.52 0.04 0.775 -0.18 -0.29
SWLS.3_recoded 0.66 0.52 0.14 0.095 . -0.28 -0.40
SWLS.4_recoded 0.60 0.52 0.08 0.441 -0.04 -0.16
SWLS.5_recoded 0.52 0.54 0.02 0.775 0.12 0.00
Code
#re-analyzing loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(SWLS_data_rasch_recoded_na_omit)

Code
#re-analysis of residual correlations
resid_cor_SWLS_recoded <- RIgetResidCor(SWLS_data_rasch_recoded_na_omit,
                                        iterations = 800,
                                        cpu = 10)

RIresidcorr(SWLS_data_rasch_recoded_na_omit, cutoff = resid_cor_SWLS_recoded$p99) 
SWLS.1_recoded SWLS.2_recoded SWLS.3_recoded SWLS.4_recoded SWLS.5_recoded
SWLS.1_recoded
SWLS.2_recoded -0.25
SWLS.3_recoded -0.17 -0.26
SWLS.4_recoded -0.2 -0.13 0.09
SWLS.5_recoded -0.3 -0.21 -0.18 -0.36
Note:
Relative cut-off value is 0.088, which is 0.285 above the average correlation (-0.198).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(SWLS_data_rasch_recoded_na_omit) #shows only significant results
[1] "No statistically significant local dependency found."
Code
partgam_LD(as.data.frame(SWLS_data_rasch_recoded_na_omit)) #shows all output, also non-significant results
            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  SWLS.1_recoded SWLS.2_recoded -0.1184 0.1573 0.4518       1     -0.4267
2  SWLS.1_recoded SWLS.3_recoded -0.0855 0.1608 0.5949       1     -0.4006
3  SWLS.1_recoded SWLS.4_recoded -0.1408 0.1350 0.2970       1     -0.4054
4  SWLS.1_recoded SWLS.5_recoded -0.1090 0.1435 0.4474       1     -0.3903
5  SWLS.2_recoded SWLS.3_recoded -0.2996 0.1735 0.0842       1     -0.6396
6  SWLS.2_recoded SWLS.4_recoded -0.0244 0.1595 0.8785       1     -0.3371
7  SWLS.2_recoded SWLS.5_recoded  0.0389 0.1666 0.8155       1     -0.2877
8  SWLS.3_recoded SWLS.4_recoded  0.3278 0.1745 0.0603       1     -0.0142
9  SWLS.3_recoded SWLS.5_recoded  0.2555 0.1764 0.1475       1     -0.0902
10 SWLS.4_recoded SWLS.5_recoded -0.1357 0.1568 0.3867       1     -0.4430
    upper
1  0.1900
2  0.2296
3  0.1238
4  0.1722
5  0.0405
6  0.2883
7  0.3654
8  0.6698
9  0.6012
10 0.1716

            Item1          Item2   gamma     se pvalue padj.BH sig   lower
1  SWLS.2_recoded SWLS.1_recoded -0.0817 0.1531 0.5935       1     -0.3818
2  SWLS.3_recoded SWLS.1_recoded  0.1312 0.1626 0.4196       1     -0.1874
3  SWLS.3_recoded SWLS.2_recoded -0.1429 0.1962 0.4665       1     -0.5274
4  SWLS.4_recoded SWLS.1_recoded  0.0678 0.1708 0.6914       1     -0.2669
5  SWLS.4_recoded SWLS.2_recoded  0.0311 0.1788 0.8618       1     -0.3193
6  SWLS.4_recoded SWLS.3_recoded  0.1508 0.2009 0.4529       1     -0.2429
7  SWLS.5_recoded SWLS.1_recoded -0.0347 0.1392 0.8029       1     -0.3077
8  SWLS.5_recoded SWLS.2_recoded  0.1048 0.1572 0.5047       1     -0.2032
9  SWLS.5_recoded SWLS.3_recoded  0.1111 0.1775 0.5313       1     -0.2367
10 SWLS.5_recoded SWLS.4_recoded -0.1698 0.1484 0.2524       1     -0.4606
    upper
1  0.2183
2  0.4499
3  0.2417
4  0.4025
5  0.3815
6  0.5445
7  0.2382
8  0.4129
9  0.4590
10 0.1210
Code
#re-analysis of targeting
RItargeting(SWLS_data_rasch_recoded_na_omit)

Code
#re-analysis of item hierarchy
itemlabels <- data.frame(itemnr = names(SWLS_data_rasch_recoded_na_omit), item = NA)

RIitemHierarchy(SWLS_data_rasch_recoded_na_omit) 

Code
#threshold locations 
RIitemparams(SWLS_data_rasch_recoded_na_omit) 
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Item location
SWLS.1_recoded -2.10 -0.53 0.64 3.02 0.25
SWLS.2_recoded -1.84 0.15 1.16 NA -0.18
SWLS.3_recoded -1.69 -0.22 1.07 NA -0.28
SWLS.4_recoded -1.14 -0.13 1.15 NA -0.04
SWLS.5_recoded -1.26 -0.23 0.17 1.80 0.12
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(SWLS_data_rasch_recoded_na_omit))

Code
RItif(na.omit(SWLS_data_rasch_recoded_na_omit), samplePSI = TRUE)

SWLS.3 was a candidate for removal but after adjusting response categories it performed slightly better (at least in the gamma coefficient analysis), so it is kept in the scale for the validity analysis.

4.3 SWLS and BEA validity analysis

Validity between SWLS and BEA is investigated by correlation analysis. As both scales has undergone Rasch analysis, participants theta estimates (computed with weighted likelihood estimation - WLE) will be correlated, and since thetas are interval data Pearson’s correlation is applied.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item data frames
BE_validity_analyses <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

SWLS_validity_analyses <- data_complete_cases[, c(
  "SWLS.1_recoded",
  "SWLS.2_recoded",
  "SWLS.3_recoded",
  "SWLS.4_recoded",
  "SWLS.5_recoded"
)]

#extracting thetas
BE_thetas <- RIestThetas(BE_validity_analyses)

SWLS_thetas <- RIestThetas(SWLS_validity_analyses)

#adding thetas to the data frame
data_complete_cases$BE_thetas_WLE <- BE_thetas$WLE 
data_complete_cases$SWLS_thetas_WLE <- SWLS_thetas$WLE
Code
correlation_BE_SWLS <- cor_test(
  data_complete_cases,
  x = "BE_thetas_WLE",
  y = "SWLS_thetas_WLE",
  method = "pearson"
)

print(correlation_BE_SWLS)
Parameter1    |      Parameter2 |    r |       95% CI | t(112) |         p
--------------------------------------------------------------------------
BE_thetas_WLE | SWLS_thetas_WLE | 0.57 | [0.43, 0.68] |   7.25 | < .001***

Observations: 114
Code
plot(correlation_BE_SWLS)

Some participants did not complete BEA and SWLS during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#data have previously been prepared for the validity sensitivity analysis
#looking only at participants who completed Bulls Eye and the other scales during the same day (N = 99)
#BE items and BE thetas for these particpants have also been prepared previously

#sensitivity analysis SWLS and BE
BE_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

SWLS_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "SWLS.1_recoded",
  "SWLS.2_recoded",
  "SWLS.3_recoded",
  "SWLS.4_recoded",
  "SWLS.5_recoded"
)]

#extracting thetas
BE_thetas_less1day <- RIestThetas(BE_validity_analyses_less1day)

SWLS_thetas_less1day <- RIestThetas(SWLS_validity_analyses_less1day)

#adding thetas to the data frame
data_complete_cases_less_than_1day$BE_thetas_less1day_WLE <- BE_thetas_less1day$WLE
data_complete_cases_less_than_1day$SWLS_thetas_less1day_WLE <- SWLS_thetas_less1day$WLE

#analysis only with participants completing Bulls Eye and ADAQ the same day
correlation_BE_SWLS_less1day <- cor_test(
  data_complete_cases_less_than_1day,
  x = "BE_thetas_less1day_WLE",
  y = "SWLS_thetas_less1day_WLE",
  method = "pearson"
)

print(correlation_BE_SWLS_less1day)
Parameter1             |               Parameter2 |    r |       95% CI | t(97) |         p
-------------------------------------------------------------------------------------------
BE_thetas_less1day_WLE | SWLS_thetas_less1day_WLE | 0.54 | [0.39, 0.67] |  6.40 | < .001***

Observations: 99
Code
#we can also check if there is a statistically significant difference 
#by adding same or different day as a predictor 
data_complete_cases_sensitivity_analysis <- data_complete_cases %>% 
  mutate(group = factor(ifelse(time_latency_survey1_bulls_eye1_days < 1,"same_day","different_day")))

m0_swls <- lm(SWLS_thetas_WLE ~ BE_thetas_WLE, data = data_complete_cases_sensitivity_analysis)
m1_swls <- lm(SWLS_thetas_WLE ~ BE_thetas_WLE + group, data = data_complete_cases_sensitivity_analysis)

anova(m0_swls, m1_swls)
Analysis of Variance Table

Model 1: SWLS_thetas_WLE ~ BE_thetas_WLE
Model 2: SWLS_thetas_WLE ~ BE_thetas_WLE + group
  Res.Df   RSS Df  Sum of Sq     F Pr(>F)
1    112 138.2                           
2    111 138.2  1 0.00078697 6e-04   0.98

Adding whether BE and SWLS were completed during the same day or not as a predictor did not significantly alter the results.

5 Valuing Questionnaire (VQ): Rasch analysis and BEA validity investigation

The Valuing Questionnaire (VQ) instrument is used as part of the validity investigation of BEA. The original VQ study (Smout et al., 2014) states that the scale has two dimensions: (A) progress in valued living (items 3, 4, 5, 7, 9), and (B) obstruction to valued living (items 1, 2, 6, 8, 10). VQ_progress and VQ_obstruction are therefore analyzed separate with Rasch analysis. Validity analyses by correlations of thetas (BEA and VQ_progress, and BEA and VQ_obstruction respectively) will be undertaken. Note that VQ was not completed by the participants from the intervention study (Hockeystudien), therefore more missing values than in the BEA Rasch analysis are to be expected.

#VQ_progess: Rasch analysis and BEA validity investigation

Code
#preparing data frame with VQ_progress items
VQ_progress_data_rasch <- bulls_eye_merged_superfinal_data[, c("VQ.3", "VQ.4", "VQ.5", "VQ.7", "VQ.9")]

#investigate missing values
RImissing(VQ_progress_data_rasch)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete VQ
VQ_progress_data_rasch_na_omit <- VQ_progress_data_rasch %>%
  select(VQ.3,
         VQ.4,
         VQ.5,
         VQ.7,
         VQ.9) %>%
  na.omit()

RImissing(VQ_progress_data_rasch_na_omit)
[1] "No missing data."
Code
#tile plot (response category frequency for the items)
RItileplot(VQ_progress_data_rasch_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(VQ_progress_data_rasch_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue  sig 
overall 61.4 29 0.00042  ***
Code
#conditional item fit (assessing unidimensionaliity)
simfit1_VQ_progress <- RIgetfit(VQ_progress_data_rasch_na_omit,
                       iterations = 100,
                       cpu = 10) 

RIitemfit(VQ_progress_data_rasch_na_omit, simfit1_VQ_progress, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
VQ.3 2.058 [0.708, 1.237] 2.519 [0.746, 1.279] 0.821 1.24 -0.26
VQ.4 0.818 [0.792, 1.174] 0.757 [0.771, 1.241] no misfit 0.014 -0.72
VQ.5 0.731 [0.793, 1.208] 0.764 [0.797, 1.238] 0.062 0.033 -0.37
VQ.7 0.698 [0.86, 1.221] 0.72 [0.85, 1.235] 0.162 0.13 -0.75
VQ.9 0.777 [0.785, 1.176] 0.812 [0.773, 1.23] 0.008 no misfit -0.71
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 71 simulated datasets.
Code
#conditional item characteristic curves
ICCplot(
  as.data.frame(VQ_progress_data_rasch_na_omit),
  itemnumber = 1:4,
  itemdescrip = c("VQ.3", "VQ.4", "VQ.5", "VQ.7"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(VQ_progress_data_rasch_na_omit),
  itemnumber = 5,
  itemdescrip = c("VQ.9"),
  method = "cut"
)

Code
#PCA of residuals (assessing unidimensionality)
RIpcmPCA(VQ_progress_data_rasch_na_omit) 
Eigenvalues Proportion of variance
1.87 55%
1.26 18.2%
1.08 15.6%
0.78 11%
0.01 0.2%
Code
#item restscores
RIrestscore(VQ_progress_data_rasch_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
VQ.3 0.12 0.45 0.33 0.000 *** 0.30 -0.26
VQ.4 0.56 0.44 0.12 0.075 . -0.15 -0.72
VQ.5 0.58 0.44 0.14 0.033 * 0.19 -0.37
VQ.7 0.61 0.42 0.19 0.001 *** -0.19 -0.75
VQ.9 0.56 0.43 0.13 0.070 . -0.15 -0.71
Code
#loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(VQ_progress_data_rasch_na_omit)

Code
#residual correlations
resid_cor_VQ_progress <- RIgetResidCor(VQ_progress_data_rasch_na_omit,
                                       iterations = 800,
                                       cpu = 10)

RIresidcorr(VQ_progress_data_rasch_na_omit, cutoff = resid_cor_VQ_progress$p99)
VQ.3 VQ.4 VQ.5 VQ.7 VQ.9
VQ.3
VQ.4 -0.6
VQ.5 -0.35 -0.07
VQ.7 -0.4 0.22 0.1
VQ.9 -0.41 0.14 -0.02 -0.08
Note:
Relative cut-off value is 0.158, which is 0.305 above the average correlation (-0.146).
Correlations above the cut-off are highlighted in red text.
Code
#partial gamma coefficients (assessing local independence)
RIpartgamLD(VQ_progress_data_rasch_na_omit) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
VQ.7 VQ.4 0.391 0.128 0.141 0.642 0.044
VQ.7 VQ.5 0.368 0.108 0.156 0.580 0.013
Code
partgam_LD(as.data.frame(VQ_progress_data_rasch_na_omit)) #shows all output, also non-significant results
   Item1 Item2   gamma     se pvalue padj.BH  sig   lower   upper
1   VQ.3  VQ.4 -0.5596 0.1072 0.0000  0.0000  *** -0.7697 -0.3494
2   VQ.3  VQ.5 -0.2258 0.1125 0.0448  0.8965      -0.4464 -0.0052
3   VQ.3  VQ.7 -0.4400 0.1284 0.0006  0.0122    * -0.6916 -0.1884
4   VQ.3  VQ.9 -0.3835 0.1252 0.0022  0.0437    * -0.6288 -0.1382
5   VQ.4  VQ.5  0.1811 0.1387 0.1915  1.0000      -0.0907  0.4529
6   VQ.4  VQ.7  0.4236 0.1443 0.0033  0.0668    .  0.1407  0.7065
7   VQ.4  VQ.9  0.3306 0.1311 0.0116  0.2329       0.0738  0.5875
8   VQ.5  VQ.7  0.3054 0.1415 0.0309  0.6188       0.0280  0.5829
9   VQ.5  VQ.9  0.1368 0.1366 0.3169  1.0000      -0.1310  0.4045
10  VQ.7  VQ.9  0.1598 0.1649 0.3323  1.0000      -0.1633  0.4829

   Item1 Item2   gamma     se pvalue padj.BH sig   lower   upper
1   VQ.4  VQ.3 -0.3189 0.1342 0.0174  0.3490     -0.5819 -0.0560
2   VQ.5  VQ.3  0.2364 0.1216 0.0519  1.0000     -0.0020  0.4747
3   VQ.5  VQ.4  0.1970 0.1374 0.1516  1.0000     -0.0722  0.4662
4   VQ.7  VQ.3  0.0000 0.1503 1.0000  1.0000     -0.2947  0.2947
5   VQ.7  VQ.4  0.3913 0.1278 0.0022  0.0439   *  0.1409  0.6417
6   VQ.7  VQ.5  0.3684 0.1081 0.0007  0.0132   *  0.1565  0.5804
7   VQ.9  VQ.3  0.0294 0.1459 0.8402  1.0000     -0.2565  0.3153
8   VQ.9  VQ.4  0.3252 0.1296 0.0121  0.2426      0.0711  0.5793
9   VQ.9  VQ.5  0.0315 0.1443 0.8272  1.0000     -0.2513  0.3143
10  VQ.9  VQ.7  0.0171 0.1635 0.9167  1.0000     -0.3033  0.3375
Code
#targeting
RItargeting(VQ_progress_data_rasch_na_omit)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(VQ_progress_data_rasch_na_omit), item = NA)

RIitemHierarchy(VQ_progress_data_rasch_na_omit) 

Code
#threshold locations 
RIitemparams(VQ_progress_data_rasch_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Threshold 6 Item location
VQ.3 0.23 -0.18 -0.58 -0.03 1.24 1.12 0.3
VQ.4 -1.16 -1.32 0.17 0.06 0.71 0.62 -0.15
VQ.5 -0.64 -0.44 -0.48 0.48 0.81 1.42 0.19
VQ.7 -1.70 -1.10 -0.50 -0.18 0.86 1.48 -0.19
VQ.9 -1.18 -1.08 -0.15 -0.01 0.70 0.84 -0.15
Note:
Item location is the average of the thresholds for each item.
Code
#analysis of response categories
RIitemCats(VQ_progress_data_rasch_na_omit, items = "all", xlims = c(-5, 5), legend = "left") 

Code
#analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(VQ_progress_data_rasch_na_omit))

Code
RItif(na.omit(VQ_progress_data_rasch_na_omit), samplePSI = TRUE)

5.1 Re-analysis of VQ_progress recoded version

Based on the analysis of the response categories and threshold locations, scale steps were merged and recoded so that each kept scale step had the highest response probability at some point on the person location scale, and also so that all threshold locations were in the right order for all items. Recoded items were given new names with the abbreviation “_recoded”.

Code
#add recoded VQ_progress items

bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(
    # Recoding VQ.3 (0+1; 2+3)
    VQ.3_recoded = recode(VQ.3,
                            `0` = 0, `1` = 0, 
                            `2` = 1, `3` = 1, 
                            `4` = 2,
                            `5` = 3,
                            `6` = 4),
    # Recoding VQ.4 (0+1; 3+4; 5+6)
    VQ.4_recoded = recode(VQ.4,
                            `0` = 0, `1` = 0, 
                            `2` = 1, 
                            `3` = 2, `4` = 2,
                            `5` = 3, `6` = 3),
    # Recoding VQ.5 (0+1+2)
    VQ.5_recoded = recode(VQ.5,
                            `0` = 0, `1` = 0,`2` = 0,
                            `3` = 1,
                            `4` = 2,
                            `5` = 3,
                            `6` = 4),
    # Recoding VQ.7 (0+1)
    VQ.7_recoded = recode(VQ.7,
                           `0` = 0, `1` = 0, 
                           `2` = 1,
                           `3` = 2, 
                           `4` = 3,
                           `5` = 4,
                           `6` = 5),
    # Recoding VQ.9 (0+1)
    VQ.9_recoded = recode(VQ.9,
                           `0` = 0, `1` = 0, 
                           `2` = 1, 
                           `3` = 2, 
                           `4` = 3,
                           `5` = 4, 
                           `6` = 5))
Code
#preparing items for VQ_progress recoded
VQ_progress_data_rasch_recoded <- bulls_eye_merged_superfinal_data[, c(
  "VQ.3_recoded",
  "VQ.4_recoded",
  "VQ.5_recoded",
  "VQ.7_recoded",
  "VQ.9_recoded"
)]

#investigate missing values
RImissing(VQ_progress_data_rasch_recoded)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete VQ
VQ_progress_data_rasch_recoded_na_omit <- VQ_progress_data_rasch_recoded %>%
  select(
    VQ.3_recoded,
    VQ.4_recoded,
    VQ.5_recoded,
    VQ.7_recoded,
    VQ.9_recoded
  ) %>%
  na.omit()

RImissing(VQ_progress_data_rasch_recoded_na_omit)
[1] "No missing data."
Code
#re-analysis of response categories using VQ_progress recoded
RIitemCats(
  VQ_progress_data_rasch_recoded_na_omit,
  items = "all",
  xlims = c(-5, 5),
  legend = "left"
)

Code
#tile plot (response category frequency for the items)
RItileplot(VQ_progress_data_rasch_recoded_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(VQ_progress_data_rasch_recoded_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 40.8 20 0.004   **
Code
#re-analysis of conditional item fit (assessing unidimensionality)
simfit1_VQ_progress_recoded <- RIgetfit(VQ_progress_data_rasch_recoded_na_omit,
                               iterations = 100,
                               cpu = 10)

RIitemfit(VQ_progress_data_rasch_recoded_na_omit, simfit1_VQ_progress_recoded, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
VQ.3_recoded 1.821 [0.826, 1.273] 1.898 [0.807, 1.282] 0.548 0.616 -0.02
VQ.4_recoded 0.721 [0.78, 1.223] 0.657 [0.751, 1.259] 0.059 0.094 -1.00
VQ.5_recoded 0.804 [0.76, 1.266] 0.804 [0.734, 1.271] no misfit no misfit 0.15
VQ.7_recoded 0.762 [0.747, 1.257] 0.768 [0.74, 1.27] no misfit no misfit -0.47
VQ.9_recoded 0.947 [0.758, 1.262] 0.944 [0.777, 1.296] no misfit no misfit -0.50
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(VQ_progress_data_rasch_recoded_na_omit),
  itemnumber = 1:4,
  itemdescrip = c(
    "VQ.3_recoded",
    "VQ.4_recoded",
    "VQ.5_recoded",
    "VQ.7_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(VQ_progress_data_rasch_recoded_na_omit),
  itemnumber = 5,
  itemdescrip = c("VQ.9_recoded"),
  method = "cut"
)

Code
#re-analysis of PCA of residuals (assessing unidimensionality)
RIpcmPCA(VQ_progress_data_rasch_recoded_na_omit)
Eigenvalues Proportion of variance
1.74 48.3%
1.37 22.6%
1.11 17.5%
0.76 11.2%
0.01 0.3%
Code
#re-analysis of item restscores
RIrestscore(VQ_progress_data_rasch_recoded_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
VQ.3_recoded 0.13 0.48 0.35 0.000 *** 0.32 -0.02
VQ.4_recoded 0.66 0.43 0.23 0.000 *** -0.65 -1.00
VQ.5_recoded 0.61 0.48 0.13 0.043 * 0.50 0.15
VQ.7_recoded 0.61 0.48 0.13 0.023 * -0.12 -0.47
VQ.9_recoded 0.55 0.49 0.06 0.404 -0.16 -0.50
Code
#re-analysis of loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(VQ_progress_data_rasch_recoded_na_omit)

Code
#re-analysis of residual correlations
resid_cor_VQ_progress_recoded <- RIgetResidCor(VQ_progress_data_rasch_recoded_na_omit,
                                               iterations = 800,
                                               cpu = 10)

RIresidcorr(VQ_progress_data_rasch_recoded_na_omit, cutoff = resid_cor_VQ_progress_recoded$p99)
VQ.3_recoded VQ.4_recoded VQ.5_recoded VQ.7_recoded VQ.9_recoded
VQ.3_recoded
VQ.4_recoded -0.46
VQ.5_recoded -0.3 -0.06
VQ.7_recoded -0.39 0.2 0.04
VQ.9_recoded -0.41 0.13 -0.19 -0.2
Note:
Relative cut-off value is 0.165, which is 0.329 above the average correlation (-0.164).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(VQ_progress_data_rasch_recoded_na_omit) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
VQ.4_recoded VQ.7_recoded 0.566 0.126 0.320 0.812 0.000
VQ.4_recoded VQ.9_recoded 0.431 0.130 0.176 0.686 0.018
Code
partgam_LD(as.data.frame(VQ_progress_data_rasch_recoded_na_omit)) #shows all output, also non-significant results
          Item1        Item2   gamma     se pvalue padj.BH  sig   lower   upper
1  VQ.3_recoded VQ.4_recoded -0.4865 0.1267 0.0001  0.0025   ** -0.7348 -0.2382
2  VQ.3_recoded VQ.5_recoded -0.2375 0.1233 0.0540  1.0000      -0.4791  0.0041
3  VQ.3_recoded VQ.7_recoded -0.3704 0.1367 0.0067  0.1347      -0.6383 -0.1025
4  VQ.3_recoded VQ.9_recoded -0.2781 0.1264 0.0278  0.5568      -0.5259 -0.0303
5  VQ.4_recoded VQ.5_recoded  0.2165 0.1774 0.2224  1.0000      -0.1313  0.5643
6  VQ.4_recoded VQ.7_recoded  0.5660 0.1256 0.0000  0.0001  ***  0.3198  0.8123
7  VQ.4_recoded VQ.9_recoded  0.4309 0.1299 0.0009  0.0183    *  0.1762  0.6856
8  VQ.5_recoded VQ.7_recoded  0.2285 0.1519 0.1327  1.0000      -0.0693  0.5263
9  VQ.5_recoded VQ.9_recoded  0.0186 0.1430 0.8966  1.0000      -0.2617  0.2989
10 VQ.7_recoded VQ.9_recoded  0.1073 0.1498 0.4740  1.0000      -0.1864  0.4009

          Item1        Item2   gamma     se pvalue padj.BH sig   lower  upper
1  VQ.4_recoded VQ.3_recoded -0.2000 0.1486 0.1783  1.0000     -0.4912 0.0912
2  VQ.5_recoded VQ.3_recoded  0.2000 0.1286 0.1198  1.0000     -0.0520 0.4520
3  VQ.5_recoded VQ.4_recoded  0.0211 0.1989 0.9157  1.0000     -0.3689 0.4110
4  VQ.7_recoded VQ.3_recoded -0.0246 0.1553 0.8742  1.0000     -0.3290 0.2798
5  VQ.7_recoded VQ.4_recoded  0.3679 0.1463 0.0119  0.2384      0.0811 0.6547
6  VQ.7_recoded VQ.5_recoded  0.3498 0.1199 0.0035  0.0705   .  0.1148 0.5848
7  VQ.9_recoded VQ.3_recoded -0.0312 0.1572 0.8424  1.0000     -0.3393 0.2768
8  VQ.9_recoded VQ.4_recoded  0.3237 0.1368 0.0179  0.3588      0.0556 0.5917
9  VQ.9_recoded VQ.5_recoded -0.0758 0.1571 0.6296  1.0000     -0.3836 0.2321
10 VQ.9_recoded VQ.7_recoded  0.0108 0.1554 0.9448  1.0000     -0.2939 0.3154
Code
#re-analysis of targeting
RItargeting(VQ_progress_data_rasch_recoded_na_omit)

Code
#re-analysis of item hierarchy
itemlabels <- data.frame(itemnr = names(VQ_progress_data_rasch_recoded_na_omit),
                         item = NA)

RIitemHierarchy(VQ_progress_data_rasch_recoded_na_omit)

Code
#re-analysis of threshold locations 
RIitemparams(VQ_progress_data_rasch_recoded_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
VQ.3_recoded -0.78 -0.02 1.00 1.10 NA 0.32
VQ.4_recoded -1.32 -0.86 0.24 NA NA -0.65
VQ.5_recoded -0.19 0.17 0.60 1.42 NA 0.5
VQ.7_recoded -1.35 -0.88 -0.47 0.65 1.47 -0.12
VQ.9_recoded -1.21 -0.50 -0.30 0.47 0.76 -0.16
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(VQ_progress_data_rasch_recoded_na_omit))

Code
RItif(na.omit(VQ_progress_data_rasch_recoded_na_omit), samplePSI = TRUE) 

  • VQ_progress recoded resulted in ordered response categories.
  • item 7 performed slightly better and is kept.
  • item 3 performs worst and bad enough to be removed. It is also the only item with indications of not belonging to the dimension.
  • item 4 also performs bad but is kept until version without item 3 is analyzed.
  • performing Rasch re-analysis with a trimmed recoded VQ_progress without item 3.

5.2 Re-analysis of VQ_progress recoded and trimmed version

Code
#preparing items in VQ_progress recoded and trimmed version
VQ_progress_data_rasch_recoded_trimmed <- bulls_eye_merged_superfinal_data[, c(
  "VQ.4_recoded",
  "VQ.5_recoded",
  "VQ.7_recoded",
  "VQ.9_recoded"
)]

#investigate missing values
RImissing(VQ_progress_data_rasch_recoded_trimmed)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete VQ
VQ_progress_data_rasch_recoded_trimmed_na_omit <- VQ_progress_data_rasch_recoded_trimmed %>%
  select(
    VQ.4_recoded,
    VQ.5_recoded,
    VQ.7_recoded,
    VQ.9_recoded
  ) %>%
  na.omit()

RImissing(VQ_progress_data_rasch_recoded_trimmed_na_omit)
[1] "No missing data."
Code
#tile plot (response category frequency for the items)
RItileplot(VQ_progress_data_rasch_recoded_trimmed_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(VQ_progress_data_rasch_recoded_trimmed_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 17.8 16 0.34      
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_VQ_progress_recoded_trimmed <- RIgetfit(VQ_progress_data_rasch_recoded_trimmed_na_omit,
                                       iterations = 100,
                                       cpu = 10) 


RIitemfit(
  VQ_progress_data_rasch_recoded_trimmed_na_omit,
  simfit1_VQ_progress_recoded_trimmed,
  cutoff = c(0.005, 0.995)
)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
VQ.4_recoded 0.779 [0.74, 1.248] 0.708 [0.715, 1.299] no misfit 0.007 -1.35
VQ.5_recoded 1.152 [0.734, 1.251] 1.125 [0.726, 1.33] no misfit no misfit 0.23
VQ.7_recoded 0.986 [0.774, 1.263] 0.954 [0.769, 1.276] no misfit no misfit -0.66
VQ.9_recoded 1.247 [0.796, 1.266] 1.181 [0.776, 1.279] no misfit no misfit -0.71
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(VQ_progress_data_rasch_recoded_trimmed_na_omit),
  itemnumber = 1:4,
  itemdescrip = c(
    "VQ.4_recoded",
    "VQ.5_recoded",
    "VQ.7_recoded",
    "VQ.9_recoded"
  )
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
#re-analysis of PCA of residuals (assessing unidimensionality)
RIpcmPCA(VQ_progress_data_rasch_recoded_trimmed_na_omit)
Eigenvalues Proportion of variance
1.61 43.6%
1.41 34.8%
0.97 21.2%
0.01 0.3%
Code
#re-analysis of item restscores
RIrestscore(VQ_progress_data_rasch_recoded_trimmed_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
VQ.4_recoded 0.73 0.60 0.13 0.059 . -0.77 -1.35
VQ.5_recoded 0.61 0.63 0.02 0.701 0.82 0.23
VQ.7_recoded 0.66 0.63 0.03 0.670 -0.07 -0.66
VQ.9_recoded 0.58 0.64 0.06 0.670 -0.12 -0.71
Code
#re-analysis of loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(VQ_progress_data_rasch_recoded_trimmed_na_omit)

Code
#re-analysis of residual correlations
resid_cor_VQ_progress_recoded_trimmed <- RIgetResidCor(VQ_progress_data_rasch_recoded_trimmed_na_omit,
                                              iterations = 800,
                                              cpu = 10)

RIresidcorr(VQ_progress_data_rasch_recoded_trimmed_na_omit, cutoff = resid_cor_VQ_progress_recoded_trimmed$p99)
VQ.4_recoded VQ.5_recoded VQ.7_recoded VQ.9_recoded
VQ.4_recoded
VQ.5_recoded -0.33
VQ.7_recoded -0.05 -0.19
VQ.9_recoded -0.13 -0.44 -0.52
Note:
Relative cut-off value is 0.037, which is 0.313 above the average correlation (-0.276).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(VQ_progress_data_rasch_recoded_trimmed_na_omit) #shows only significant results
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
VQ.4_recoded VQ.7_recoded 0.41 0.143 0.130 0.690 0.049
VQ.4_recoded VQ.9_recoded 0.4 0.133 0.138 0.662 0.033
Code
partgam_LD(as.data.frame(VQ_progress_data_rasch_recoded_trimmed_na_omit)) #shows all output, also non-significant results 
         Item1        Item2   gamma     se pvalue padj.BH sig   lower  upper
1 VQ.4_recoded VQ.5_recoded -0.0405 0.1781 0.8203  1.0000     -0.3895 0.3086
2 VQ.4_recoded VQ.7_recoded  0.4100 0.1428 0.0041  0.0492   *  0.1300 0.6900
3 VQ.4_recoded VQ.9_recoded  0.4000 0.1335 0.0027  0.0328   *  0.1383 0.6617
4 VQ.5_recoded VQ.7_recoded  0.0042 0.1610 0.9791  1.0000     -0.3113 0.3197
5 VQ.5_recoded VQ.9_recoded -0.2245 0.1361 0.0990  1.0000     -0.4912 0.0422
6 VQ.7_recoded VQ.9_recoded -0.0705 0.1426 0.6208  1.0000     -0.3500 0.2089

         Item1        Item2   gamma     se pvalue padj.BH sig   lower   upper
1 VQ.5_recoded VQ.4_recoded -0.3333 0.1689 0.0484  0.5809     -0.6643 -0.0023
2 VQ.7_recoded VQ.4_recoded  0.2577 0.1610 0.1095  1.0000     -0.0579  0.5733
3 VQ.7_recoded VQ.5_recoded  0.1907 0.1500 0.2037  1.0000     -0.1034  0.4847
4 VQ.9_recoded VQ.4_recoded  0.1681 0.1560 0.2812  1.0000     -0.1376  0.4737
5 VQ.9_recoded VQ.5_recoded -0.1349 0.1400 0.3349  1.0000     -0.4093  0.1394
6 VQ.9_recoded VQ.7_recoded -0.2549 0.1427 0.0740  0.8885     -0.5346  0.0248
Code
#re-analysis of targeting
RItargeting(VQ_progress_data_rasch_recoded_trimmed_na_omit)

Code
#re-analysis of item hierarchy
itemlabels <- data.frame(itemnr = names(VQ_progress_data_rasch_recoded_trimmed_na_omit),
                         item = NA)

RIitemHierarchy(VQ_progress_data_rasch_recoded_trimmed_na_omit)

Code
#re-analysis of threshold locations 
RIitemparams(VQ_progress_data_rasch_recoded_trimmed_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Item location
VQ.4_recoded -1.91 -0.98 0.59 NA NA -0.77
VQ.5_recoded -0.41 0.26 1.03 2.39 NA 0.82
VQ.7_recoded -2.08 -1.15 -0.48 0.97 2.38 -0.07
VQ.9_recoded -1.81 -0.73 -0.30 0.73 1.51 -0.12
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI
RItif(na.omit(VQ_progress_data_rasch_recoded_trimmed_na_omit))

Code
RItif(na.omit(VQ_progress_data_rasch_recoded_trimmed_na_omit), samplePSI = TRUE)

  • Item 4 performs slightly better but still not optimal. It is however kept in the scale for the validity analysis.
  • The other items works OK but TIF (reliability) is not good at all.

5.3 VQ_progress and BEA validity analysis

Validity between VQ_progress and BEA is investigated by correlation analysis. As both scales has undergone Rasch analysis, participants theta estimates (computed with weighted likelihood estimation - WLE) will be correlated, and since thetas are interval data, Pearson’s correlation is applied.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item data frames
BE_validity_analyses <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

VQ_progress_validity_analyses <- data_complete_cases[, c(
  "VQ.4_recoded",
  "VQ.5_recoded",
  "VQ.7_recoded",
  "VQ.9_recoded"
)]

#extracting thetas
BE_thetas <- RIestThetas(BE_validity_analyses)

VQ_progress_thetas <- RIestThetas(VQ_progress_validity_analyses)

#adding thetas to the data frame
data_complete_cases$BE_thetas_WLE <- BE_thetas$WLE 
data_complete_cases$VQ_progress_thetas_WLE <- VQ_progress_thetas$WLE
Code
correlation_BE_VQ_progress <- cor_test(
  data_complete_cases,
  x = "BE_thetas_WLE",
  y = "VQ_progress_thetas_WLE",
  method = "pearson"
)

print(correlation_BE_VQ_progress)
Parameter1    |             Parameter2 |    r |       95% CI | t(112) |         p
---------------------------------------------------------------------------------
BE_thetas_WLE | VQ_progress_thetas_WLE | 0.52 | [0.37, 0.64] |   6.41 | < .001***

Observations: 114
Code
plot(correlation_BE_VQ_progress)

Some participants did not complete BEA and VQ during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#data have previously been prepared for the validity sensitivity analysis
#looking only at participants who completed Bulls Eye and the other scales during the same day (N = 99)
#BE items and BE thetas for these particpants have also been prepared previously

#sensitivity analysis VQ_progress and BE
BE_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

VQ_progress_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "VQ.4_recoded",
  "VQ.5_recoded",
  "VQ.7_recoded",
  "VQ.9_recoded"
)]

#extracting thetas
BE_thetas_less1day <- RIestThetas(BE_validity_analyses_less1day)

VQ_progress_thetas_less1day <- RIestThetas(VQ_progress_validity_analyses_less1day)

#adding thetas to the data frame
data_complete_cases_less_than_1day$BE_thetas_less1day_WLE <- BE_thetas_less1day$WLE

data_complete_cases_less_than_1day$VQ_progress_thetas_less1day_WLE <- VQ_progress_thetas_less1day$WLE

#analysis only with participants completing Bulls Eye and ADAQ the same day
correlation_BE_VQ_progress_less1day <- cor_test(
  data_complete_cases_less_than_1day,
  x = "BE_thetas_less1day_WLE",
  y = "VQ_progress_thetas_less1day_WLE",
  method = "pearson"
)

print(correlation_BE_VQ_progress_less1day)
Parameter1             |                      Parameter2 |    r |       95% CI
------------------------------------------------------------------------------
BE_thetas_less1day_WLE | VQ_progress_thetas_less1day_WLE | 0.51 | [0.35, 0.65]

Parameter1             | t(97) |         p
------------------------------------------
BE_thetas_less1day_WLE |  5.91 | < .001***

Observations: 99
Code
#we can also check if there is a statistically significant difference 
#by adding same or different day as a predictor 
data_complete_cases_sensitivity_analysis <- data_complete_cases %>% 
  mutate(group = factor(ifelse(time_latency_survey1_bulls_eye1_days < 1,"same_day","different_day")))

m0_VQ_progress <- lm(VQ_progress_thetas_WLE ~ BE_thetas_WLE, data = data_complete_cases_sensitivity_analysis)

m1_VQ_progress <- lm(VQ_progress_thetas_WLE ~ BE_thetas_WLE + group, data = data_complete_cases_sensitivity_analysis)

anova(m0_VQ_progress, m1_VQ_progress)
Analysis of Variance Table

Model 1: VQ_progress_thetas_WLE ~ BE_thetas_WLE
Model 2: VQ_progress_thetas_WLE ~ BE_thetas_WLE + group
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    112 202.59                              
2    111 196.38  1    6.2105 3.5103 0.06362 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding whether BE and VQ_progress were completed during the same day or not as a predictor did not significantly alter the results.

Now moving on to analyzing the other part of VQ, VQ_obstruction (to valued living, items 1, 2, 6, 8, 10).

#VQ_obstruction: Rasch analysis and BEA validity investigation

Code
#preparing data frame with VQ_obstruction items
VQ_obstruction_data_rasch <- bulls_eye_merged_superfinal_data[, c("VQ.1", "VQ.2", "VQ.6", "VQ.8", "VQ.10")]

#investigating missing values
RImissing(VQ_obstruction_data_rasch)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete VQ

VQ_obstruction_data_rasch_na_omit <- VQ_obstruction_data_rasch %>%
  select(VQ.1, VQ.2, VQ.6, VQ.8, VQ.10) %>%
  na.omit()

RImissing(VQ_obstruction_data_rasch_na_omit)
[1] "No missing data."
Code
#tile plot (response category frequency for the items)
RItileplot(VQ_obstruction_data_rasch_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(VQ_obstruction_data_rasch_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 32.9 29 0.28      
Code
#conditional item fit (assessing unidimensionality)
simfit1_VQ_obstruction <- RIgetfit(VQ_obstruction_data_rasch_na_omit,
                                   iterations = 100,
                                   cpu = 10) 

RIitemfit(VQ_obstruction_data_rasch_na_omit, simfit1_VQ_obstruction, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
VQ.1 0.927 [0.799, 1.277] 0.951 [0.811, 1.277] no misfit no misfit -0.28
VQ.2 0.981 [0.748, 1.242] 0.971 [0.744, 1.297] no misfit no misfit -0.42
VQ.6 0.989 [0.77, 1.204] 0.977 [0.741, 1.241] no misfit no misfit -0.53
VQ.8 1.038 [0.832, 1.232] 0.994 [0.784, 1.242] no misfit no misfit -1.10
VQ.10 1.091 [0.731, 1.21] 1.158 [0.745, 1.376] no misfit no misfit -0.51
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 54 simulated datasets.
Code
#conditional item characteristic curves
ICCplot(
  as.data.frame(VQ_obstruction_data_rasch_na_omit),
  itemnumber = 1:4,
  itemdescrip = c("VQ.1", "VQ.2", "VQ.6", "VQ.8"),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(VQ_obstruction_data_rasch_na_omit),
  itemnumber = 5,
  itemdescrip = c("VQ.10"),
  method = "cut"
)

Code
#PCA of residuals (assessing unidimensionality)
RIpcmPCA(VQ_obstruction_data_rasch_na_omit) 
Eigenvalues Proportion of variance
1.69 34.8%
1.27 24.4%
1.07 20.9%
0.97 19.8%
0.01 0.2%
Code
#item restscores
RIrestscore(VQ_obstruction_data_rasch_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
VQ.1 0.47 0.44 0.03 0.991 0.29 -0.28
VQ.2 0.47 0.44 0.03 0.991 0.14 -0.42
VQ.6 0.46 0.44 0.02 0.991 0.04 -0.53
VQ.8 0.40 0.39 0.01 0.991 -0.53 -1.10
VQ.10 0.45 0.45 0.00 0.991 0.06 -0.51
Code
#loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(VQ_obstruction_data_rasch_na_omit)

Code
#residual correlations
resid_cor_VQ_obstruction <- RIgetResidCor(VQ_obstruction_data_rasch_na_omit,
                                          iterations = 800,
                                          cpu = 10)

RIresidcorr(VQ_obstruction_data_rasch_na_omit, cutoff = resid_cor_VQ_obstruction$p99) 
VQ.1 VQ.2 VQ.6 VQ.8 VQ.10
VQ.1
VQ.2 -0.01
VQ.6 -0.06 -0.19
VQ.8 -0.28 -0.16 -0.3
VQ.10 -0.35 -0.35 -0.24 0.05
Note:
Relative cut-off value is 0.1, which is 0.289 above the average correlation (-0.189).
Correlations above the cut-off are highlighted in red text.
Code
#partial gamma coefficients (assessing local independence)
RIpartgamLD(VQ_obstruction_data_rasch_na_omit) #shows only significant results
[1] "No statistically significant local dependency found."
Code
partgam_LD(as.data.frame(VQ_obstruction_data_rasch_na_omit)) #shows all output, also non-significant results
   Item1 Item2   gamma     se pvalue padj.BH sig   lower  upper
1   VQ.1  VQ.2  0.1681 0.1483 0.2570       1     -0.1226 0.4587
2   VQ.1  VQ.6  0.1791 0.1224 0.1434       1     -0.0608 0.4190
3   VQ.1  VQ.8 -0.0851 0.1615 0.5982       1     -0.4016 0.2314
4   VQ.1 VQ.10 -0.1184 0.1418 0.4039       1     -0.3963 0.1596
5   VQ.2  VQ.6  0.0424 0.1339 0.7517       1     -0.2201 0.3048
6   VQ.2  VQ.8 -0.0122 0.1570 0.9381       1     -0.3199 0.2955
7   VQ.2 VQ.10  0.0404 0.1374 0.7689       1     -0.2289 0.3096
8   VQ.6  VQ.8 -0.1546 0.1473 0.2938       1     -0.4433 0.1341
9   VQ.6 VQ.10 -0.1811 0.1394 0.1938       1     -0.4543 0.0920
10  VQ.8 VQ.10  0.2162 0.1543 0.1612       1     -0.0863 0.5187

   Item1 Item2   gamma     se pvalue padj.BH sig   lower   upper
1   VQ.2  VQ.1  0.2544 0.1305 0.0512  1.0000     -0.0014  0.5101
2   VQ.6  VQ.1  0.2297 0.1243 0.0646  1.0000     -0.0139  0.4734
3   VQ.6  VQ.2 -0.0123 0.1336 0.9264  1.0000     -0.2743  0.2496
4   VQ.8  VQ.1 -0.3950 0.1319 0.0028  0.0552   . -0.6536 -0.1364
5   VQ.8  VQ.2 -0.0974 0.1480 0.5102  1.0000     -0.3875  0.1926
6   VQ.8  VQ.6 -0.1142 0.1493 0.4444  1.0000     -0.4067  0.1784
7  VQ.10  VQ.1 -0.2483 0.1308 0.0578  1.0000     -0.5047  0.0082
8  VQ.10  VQ.2 -0.0386 0.1484 0.7946  1.0000     -0.3294  0.2521
9  VQ.10  VQ.6 -0.1613 0.1430 0.2595  1.0000     -0.4416  0.1191
10 VQ.10  VQ.8  0.1919 0.1441 0.1828  1.0000     -0.0904  0.4743
Code
#targeting
RItargeting(VQ_obstruction_data_rasch_na_omit)

Code
#item hierarchy
itemlabels <- data.frame(itemnr = names(VQ_obstruction_data_rasch_na_omit),
                         item = NA)

RIitemHierarchy(VQ_obstruction_data_rasch_na_omit)

Code
#threshold locations 
RIitemparams(VQ_obstruction_data_rasch_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Threshold 6 Item location
VQ.1 -0.06 -0.75 -0.16 0.74 0.86 1.10 0.29
VQ.2 0.41 -1.85 0.05 0.58 1.24 0.42 0.14
VQ.6 -1.30 -0.60 0.39 0.54 0.20 1.01 0.04
VQ.8 -0.38 -2.13 0.02 -0.87 0.13 0.05 -0.53
VQ.10 -1.10 -0.28 0.06 0.70 0.24 0.74 0.06
Note:
Item location is the average of the thresholds for each item.
Code
#analysis of response categories
RIitemCats(
  VQ_obstruction_data_rasch_na_omit,
  items = "all",
  xlims = c(-5, 5),
  legend = "left"
)

Code
#analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(VQ_obstruction_data_rasch_na_omit))

Code
RItif(na.omit(VQ_obstruction_data_rasch_na_omit), samplePSI = TRUE)

5.4 Re-analysis of VQ_obstruction recoded version

Based on the analysis of the response categories and threshold locations, scale steps were merged and recoded so that each kept scale step had the highest response probability at some point on the person location scale, and also so that all threshold locations were in the right order for all items. Recoded items were given new names with the abbreviation “_recoded”.

Code
#add recoded VQ_obstruction items

bulls_eye_merged_superfinal_data <- bulls_eye_merged_superfinal_data %>%
  mutate(
    # Recoding VQ.1 (0+1+2)
    VQ.1_recoded = recode(VQ.1,
                            `0` = 0, `1` = 0, `2` = 0, 
                            `3` = 1,
                            `4` = 2,
                            `5` = 3,
                            `6` = 4),
    # Recoding VQ.2 (0+1; 4+5)
    VQ.2_recoded = recode(VQ.2,
                            `0` = 0, `1` = 0, 
                            `2` = 1,
                            `3` = 2, 
                            `4` = 3, `5` = 3,
                            `6` = 4),
    # Recoding VQ.6 (0+1; 3+4)
    VQ.6_recoded = recode(VQ.6,
                           `0` = 0, `1` = 0,
                           `2` = 1, 
                           `3` = 2, `4` = 2,
                           `5` = 3,
                           `6` = 4),
    # Recoding VQ.8 (0+1+2+3; 4+5)
    VQ.8_recoded = recode(VQ.8,
                           `0` = 0, `1` = 0, `2` = 0, `3` = 0, 
                           `4` = 1, `5` = 1,
                           `6` = 2),
    # Recoding VQ.10 (0+1, 2+3; 4+5)
    VQ.10_recoded = recode(VQ.10,
                           `0` = 0, `1` = 0,
                           `2` = 1, `3` = 1,
                           `4` = 2, `5` = 2,
                           `6` = 3))
Code
VQ_obstruction_data_rasch_recoded <- bulls_eye_merged_superfinal_data[, c(
  "VQ.1_recoded",
  "VQ.2_recoded",
  "VQ.6_recoded",
  "VQ.8_recoded",
  "VQ.10_recoded"
)]

#investigating missing values
RImissing(VQ_obstruction_data_rasch_recoded)

Code
#removing missing values 
#in practice this means removing participants from the intervention study (Hockeystudien) that did not complete VQ
VQ_obstruction_data_rasch_recoded_na_omit <- VQ_obstruction_data_rasch_recoded %>%
  select(
    VQ.1_recoded,
    VQ.2_recoded,
    VQ.6_recoded,
    VQ.8_recoded,
    VQ.10_recoded
  ) %>%
  na.omit()

RImissing(VQ_obstruction_data_rasch_recoded_na_omit)
[1] "No missing data."
Code
#re-analysis of response categories using VQ_obstruction recoded
RIitemCats(
  VQ_obstruction_data_rasch_recoded_na_omit,
  items = "all",
  xlims = c(-5, 5),
  legend = "left"
)

Code
#tile plot (response category frequency for the items)
RItileplot(VQ_obstruction_data_rasch_recoded_na_omit)

Code
#conditional likelihood ratio tests (CLR)
#a test of homogeneity that assesses whether item parameters are homogeneous across all levels of theta 
clr_tests(VQ_obstruction_data_rasch_recoded_na_omit, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig
overall 21.9 16 0.15      
Code
#re-analysis of conditional item fit (assessing unidimensionaliity)
simfit1_VQ_obstruction_recoded <- RIgetfit(VQ_obstruction_data_rasch_recoded_na_omit,
                                           iterations = 100,
                                           cpu = 10)

RIitemfit(VQ_obstruction_data_rasch_recoded_na_omit, simfit1_VQ_obstruction_recoded, cutoff = c(0.005, 0.995))
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Relative location
VQ.1_recoded 1.091 [0.793, 1.177] 1.01 [0.792, 1.207] no misfit no misfit 0.23
VQ.2_recoded 0.997 [0.762, 1.356] 0.995 [0.755, 1.395] no misfit no misfit -0.47
VQ.6_recoded 0.968 [0.761, 1.32] 0.945 [0.769, 1.331] no misfit no misfit -0.27
VQ.8_recoded 1.039 [0.798, 1.175] 1.003 [0.8, 1.211] no misfit no misfit -0.79
VQ.10_recoded 1.016 [0.769, 1.228] 1.05 [0.744, 1.331] no misfit no misfit -0.31
Note:
MSQ values based on conditional calculations (n = 122 complete cases).
Simulation based thresholds from 100 simulated datasets.
Code
#re-analysis of conditional item characteristic curves
ICCplot(
  as.data.frame(VQ_obstruction_data_rasch_recoded_na_omit),
  itemnumber = 1:4,
  itemdescrip = c(
    "VQ.1_recoded",
    "VQ.2_recoded",
    "VQ.6_recoded",
    "VQ.8_recoded"
  ),
  method = "cut"
)

[1] "Please press Zoom on the Plots window to see the plot"
Code
ICCplot(
  as.data.frame(VQ_obstruction_data_rasch_recoded_na_omit),
  itemnumber = 5,
  itemdescrip = c("VQ.10_recoded"),
  method = "cut"
)

Code
#re-analyzing PCA of residuals (assessing unidimensionality)
RIpcmPCA(VQ_obstruction_data_rasch_recoded_na_omit)
Eigenvalues Proportion of variance
1.61 33.4%
1.32 25.1%
1.17 22.4%
0.88 19%
0.01 0.1%
Code
#re-analysis of item restscores
RIrestscore(VQ_obstruction_data_rasch_recoded_na_omit)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
VQ.1_recoded 0.50 0.49 0.01 0.926 0.50 0.23
VQ.2_recoded 0.50 0.46 0.04 0.926 -0.20 -0.47
VQ.6_recoded 0.48 0.47 0.01 0.926 0.00 -0.27
VQ.8_recoded 0.37 0.40 0.03 0.926 -0.52 -0.79
VQ.10_recoded 0.46 0.44 0.02 0.926 -0.05 -0.31
Code
#re-analysis of loadings on the first residual contrast (assessing clusters in data or multidimensionality)
RIloadLoc(VQ_obstruction_data_rasch_recoded_na_omit)

Code
#re-analysis of residual correlations
resid_cor_VQ_obstruction_recoded <- RIgetResidCor(VQ_obstruction_data_rasch_recoded_na_omit,
                                                  iterations = 800,
                                                  cpu = 10)

RIresidcorr(VQ_obstruction_data_rasch_recoded_na_omit, cutoff = resid_cor_VQ_obstruction_recoded$p99)
VQ.1_recoded VQ.2_recoded VQ.6_recoded VQ.8_recoded VQ.10_recoded
VQ.1_recoded
VQ.2_recoded -0.15
VQ.6_recoded -0.16 -0.25
VQ.8_recoded -0.35 -0.14 -0.15
VQ.10_recoded -0.29 -0.22 -0.23 0.13
Note:
Relative cut-off value is 0.153, which is 0.335 above the average correlation (-0.182).
Correlations above the cut-off are highlighted in red text.
Code
#re-analysis of partial gamma coefficients (assessing local independence)
RIpartgamLD(VQ_obstruction_data_rasch_recoded_na_omit) #shows only significant results
[1] "No statistically significant local dependency found."
Code
partgam_LD(as.data.frame(VQ_obstruction_data_rasch_recoded_na_omit)) #shows all output, also non-significant results
          Item1         Item2   gamma     se pvalue padj.BH sig   lower  upper
1  VQ.1_recoded  VQ.2_recoded  0.1521 0.1423 0.2852       1     -0.1268 0.4311
2  VQ.1_recoded  VQ.6_recoded  0.1754 0.1309 0.1803       1     -0.0811 0.4319
3  VQ.1_recoded  VQ.8_recoded -0.2906 0.1659 0.0797       1     -0.6157 0.0345
4  VQ.1_recoded VQ.10_recoded -0.1321 0.1571 0.4005       1     -0.4400 0.1758
5  VQ.2_recoded  VQ.6_recoded -0.0168 0.1382 0.9031       1     -0.2878 0.2541
6  VQ.2_recoded  VQ.8_recoded -0.0806 0.1669 0.6293       1     -0.4077 0.2466
7  VQ.2_recoded VQ.10_recoded  0.0000 0.1464 1.0000       1     -0.2869 0.2869
8  VQ.6_recoded  VQ.8_recoded -0.0132 0.1632 0.9355       1     -0.3331 0.3066
9  VQ.6_recoded VQ.10_recoded -0.1055 0.1549 0.4960       1     -0.4090 0.1981
10 VQ.8_recoded VQ.10_recoded  0.2750 0.1681 0.1019       1     -0.0545 0.6045

           Item1        Item2   gamma     se pvalue padj.BH sig   lower   upper
1   VQ.2_recoded VQ.1_recoded  0.1180 0.1310 0.3679  1.0000     -0.1388  0.3748
2   VQ.6_recoded VQ.1_recoded  0.1599 0.1241 0.1977  1.0000     -0.0834  0.4032
3   VQ.6_recoded VQ.2_recoded -0.0789 0.1373 0.5654  1.0000     -0.3481  0.1903
4   VQ.8_recoded VQ.1_recoded -0.3433 0.1473 0.0198  0.3953     -0.6320 -0.0546
5   VQ.8_recoded VQ.2_recoded -0.0168 0.1709 0.9217  1.0000     -0.3518  0.3182
6   VQ.8_recoded VQ.6_recoded -0.0866 0.1577 0.5828  1.0000     -0.3957  0.2224
7  VQ.10_recoded VQ.1_recoded -0.1026 0.1452 0.4799  1.0000     -0.3871  0.1820
8  VQ.10_recoded VQ.2_recoded -0.0609 0.1648 0.7116  1.0000     -0.3839  0.2621
9  VQ.10_recoded VQ.6_recoded -0.1447 0.1508 0.3372  1.0000     -0.4403  0.1509
10 VQ.10_recoded VQ.8_recoded  0.3488 0.1661 0.0357  0.7136      0.0233  0.6743
Code
#re-analysis of targeting
RItargeting(VQ_obstruction_data_rasch_recoded_na_omit)

Code
#re-analysis of item hierarchy
itemlabels <- data.frame(itemnr = names(VQ_obstruction_data_rasch_recoded_na_omit),
                         item = NA)

RIitemHierarchy(VQ_obstruction_data_rasch_recoded_na_omit)

Code
#threshold locations 
RIitemparams(VQ_obstruction_data_rasch_recoded_na_omit)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Item location
VQ.1_recoded -0.08 0.38 0.68 1.02 0.5
VQ.2_recoded -1.47 -0.36 -0.11 1.13 -0.2
VQ.6_recoded -0.83 -0.60 0.58 0.83 0
VQ.8_recoded -1.29 0.25 NA NA -0.52
VQ.10_recoded -1.14 0.06 0.94 NA -0.05
Note:
Item location is the average of the thresholds for each item.
Code
#re-analysis of reliability with Test Information Function (TIF)
#first output shows without sample PSI and second output with sample PSI 
RItif(na.omit(VQ_obstruction_data_rasch_recoded_na_omit))

Code
RItif(na.omit(VQ_obstruction_data_rasch_recoded_na_omit), samplePSI = TRUE)

Response categories in VQ_obstruction is now ordered. All items in VQ_obstruction works fine enough for all to be kept. Moving on to the validity analysis.

5.5 VQ_obstruction and BEA validity analysis

Validity between VQ_obstruction and BEA is investigated by correlation analysis. As both scales has undergone Rasch analysis, participants theta estimates (computed with weighted likelihood estimation - WLE) will be correlated, and since thetas are interval data, Pearson’s correlation is applied.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item data frames
BE_validity_analyses <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

VQ_obstruction_validity_analyses <- data_complete_cases[, c(
  "VQ.1_recoded",
  "VQ.2_recoded",
  "VQ.6_recoded",
  "VQ.8_recoded",
  "VQ.10_recoded"
)]

#extracting thetas
BE_thetas <- RIestThetas(BE_validity_analyses)

VQ_obstruction_thetas <- RIestThetas(VQ_obstruction_validity_analyses)

#adding thetas to the data frame
data_complete_cases$BE_thetas_WLE <- BE_thetas$WLE 
data_complete_cases$VQ_obstruction_thetas_WLE <- VQ_obstruction_thetas$WLE
Code
correlation_BE_VQ_obstruction <- cor_test(
  data_complete_cases,
  x = "BE_thetas_WLE",
  y = "VQ_obstruction_thetas_WLE",
  method = "pearson"
)

print(correlation_BE_VQ_obstruction)
Parameter1    |                Parameter2 |    r |       95% CI | t(112) |         p
------------------------------------------------------------------------------------
BE_thetas_WLE | VQ_obstruction_thetas_WLE | 0.45 | [0.29, 0.59] |   5.36 | < .001***

Observations: 114
Code
plot(correlation_BE_VQ_obstruction)

Some participants did not complete BEA and VQ during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#data have previously been prepared for the validity sensitivity analysis
#looking only at participants who completed Bulls Eye and the other scales during the same day (N = 99)
#BE items and BE thetas for these particpants have also been prepared previously

#sensitivity analysis VQ_obstruction and BE
BE_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

VQ_obstruction_validity_analyses_less1day <- data_complete_cases_less_than_1day[, c(
  "VQ.1_recoded",
  "VQ.2_recoded",
  "VQ.6_recoded",
  "VQ.8_recoded",
  "VQ.10_recoded"
)]

#extracting thetas
BE_thetas_less1day <- RIestThetas(BE_validity_analyses_less1day)

VQ_obstruction_thetas_less1day <- RIestThetas(VQ_obstruction_validity_analyses_less1day)

#adding thetas to the data frame
data_complete_cases_less_than_1day$BE_thetas_less1day_WLE <- BE_thetas_less1day$WLE

data_complete_cases_less_than_1day$VQ_obstruction_thetas_less1day_WLE <- VQ_obstruction_thetas_less1day$WLE

#analysis only with participants completing Bulls Eye and ADAQ the same day
correlation_BE_VQ_obstruction_less1day <- cor_test(
  data_complete_cases_less_than_1day,
  x = "BE_thetas_less1day_WLE",
  y = "VQ_obstruction_thetas_less1day_WLE",
  method = "pearson"
)

print(correlation_BE_VQ_obstruction_less1day)
Parameter1             |                         Parameter2 |    r
------------------------------------------------------------------
BE_thetas_less1day_WLE | VQ_obstruction_thetas_less1day_WLE | 0.42

Parameter1             |       95% CI | t(97) |         p
---------------------------------------------------------
BE_thetas_less1day_WLE | [0.24, 0.57] |  4.51 | < .001***

Observations: 99
Code
#we can also check if there is a statistically significant difference 
#by adding same or different day as a predictor 
data_complete_cases_sensitivity_analysis <- data_complete_cases %>% 
  mutate(group = factor(ifelse(time_latency_survey1_bulls_eye1_days < 1,"same_day","different_day")))

m0_VQ_obstruction <- lm(VQ_obstruction_thetas_WLE ~ BE_thetas_WLE, data = data_complete_cases_sensitivity_analysis)

m1_VQ_obstruction <- lm(VQ_obstruction_thetas_WLE ~ BE_thetas_WLE + group, data = data_complete_cases_sensitivity_analysis)

anova(m0_VQ_obstruction, m1_VQ_obstruction)
Analysis of Variance Table

Model 1: VQ_obstruction_thetas_WLE ~ BE_thetas_WLE
Model 2: VQ_obstruction_thetas_WLE ~ BE_thetas_WLE + group
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    112 85.384                           
2    111 85.360  1  0.023474 0.0305 0.8616

Adding whether BE and VQ_obstruction were completed during the same day or not as a predictor did not significantly alter the results.

6 Sport Motivation Scale-2 and BEA validity analysis

Validity between Sport Motivation Scale-2 (SMS-2) and BEA is investigated by correlation analysis. SMS-2 has not undergone Rasch analysis, thus participants sum scores will be correlated. Since sum scores are ordinal data, Spearman’s rank correlation is applied.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item preparation and plot of sum scores 
SMS.2_items <- data_complete_cases[, c(
  "SMS-2.1",
  "SMS-2.2",
  "SMS-2.3",
  "SMS-2.4",
  "SMS-2.5",
  "SMS-2.6",
  "SMS-2.7",
  "SMS-2.8",
  "SMS-2.9",
  "SMS-2.10",
  "SMS-2.11",
  "SMS-2.12",
  "SMS-2.13",
  "SMS-2.14",
  "SMS-2.15",
  "SMS-2.16",
  "SMS-2.17",
  "SMS-2.18"
)]


BE_items <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

BE_items_names <- names(BE_items)
SMS.2_items_names <- names(SMS.2_items)

#computing sum scores for each scale
data_complete_cases_SMS.2_BE <- data_complete_cases %>%
  mutate(SMS_2_sum = rowSums(select(., all_of(SMS.2_items_names)), na.rm = TRUE),
         BE_sum = rowSums(select(., all_of(BE_items_names)), na.rm = TRUE))

#plot of sum scores (raw data, not ranked)
ggplot(data_complete_cases_SMS.2_BE, aes(x = SMS_2_sum, y = BE_sum)) +
  geom_point() +
  labs(
    title = "Scatterplot of SMS-2 vs BE Sum Scores",
    x = "SMS-2 Sum Score",
    y = "BE Sum Score"
  ) +
  theme_minimal()

Code
#compute Spearman's rank correlation
correlation_SMS.2_BE <- correlation(data_complete_cases_SMS.2_BE[, c("SMS_2_sum", "BE_sum")], method = "spearman")
print(correlation_SMS.2_BE)
# Correlation Matrix (spearman-method)

Parameter1 | Parameter2 |  rho |        95% CI |        S |     p
-----------------------------------------------------------------
SMS_2_sum  |     BE_sum | 0.11 | [-0.09, 0.29] | 2.21e+05 | 0.264

p-value adjustment method: Holm (1979)
Observations: 114

Some participants did not complete BEA and SMS-2 during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#sensitivity analysis correlation SMS-2 and BEA
SMS.2_items_less1day <- data_complete_cases_less_than_1day[, c(
  "SMS-2.1",
  "SMS-2.2",
  "SMS-2.3",
  "SMS-2.4",
  "SMS-2.5",
  "SMS-2.6",
  "SMS-2.7",
  "SMS-2.8",
  "SMS-2.9",
  "SMS-2.10",
  "SMS-2.11",
  "SMS-2.12",
  "SMS-2.13",
  "SMS-2.14",
  "SMS-2.15",
  "SMS-2.16",
  "SMS-2.17",
  "SMS-2.18"
)]


BE_items_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

BE_items_less1day_names <- names(BE_items_less1day)
SMS.2_items_less1day_names <- names(SMS.2_items_less1day)


#computing sum scores for each scale
data_complete_cases_SMS.2_BE_less1day <- data_complete_cases_less_than_1day %>%
  mutate(SMS_2_sum = rowSums(select(., all_of(
    SMS.2_items_less1day_names
  )), na.rm = TRUE),
  BE_sum = rowSums(select(., all_of(
    BE_items_less1day_names
  )), na.rm = TRUE))

#compute Spearman's rank correlation
correlation_SMS.2_BE_less1day <- correlation(data_complete_cases_SMS.2_BE_less1day[, c("SMS_2_sum", "BE_sum")], method = "spearman")

print(correlation_SMS.2_BE_less1day)
# Correlation Matrix (spearman-method)

Parameter1 | Parameter2 |  rho |        95% CI |        S |     p
-----------------------------------------------------------------
SMS_2_sum  |     BE_sum | 0.09 | [-0.12, 0.29] | 1.47e+05 | 0.384

p-value adjustment method: Holm (1979)
Observations: 99

Result from sensitivity analysis correlation is very close to the original correlation.

7 Competitive State Anxiety Inventory-2 Revised and BEA validity analysis

Validity between Competitive State Anxiety Inventory-2 Revised (CSAI-2R) and BEA is investigated by correlation analysis. CSAI-2R has not undergone Rasch analysis, thus participants sum scores will be correlated. Since sum scores are ordinal data, Spearman’s rank correlation is applied.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item preparation and plot of sum scores
CSAI2R_items <- data_complete_cases[, c(
  "CSAI-2R.1",
  "CSAI-2R.2",
  "CSAI-2R.3",
  "CSAI-2R.4",
  "CSAI-2R.5",
  "CSAI-2R.6",
  "CSAI-2R.7",
  "CSAI-2R.8",
  "CSAI-2R.9",
  "CSAI-2R.10",
  "CSAI-2R.11",
  "CSAI-2R.12",
  "CSAI-2R.13",
  "CSAI-2R.14",
  "CSAI-2R.15",
  "CSAI-2R.16",
  "CSAI-2R.17"
)]

BE_items <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

BE_items_names <- names(BE_items)
CSAI2R_items_names <- names(CSAI2R_items)

#computing sum scores for each scale
data_complete_cases_CSAI2R_BE <- data_complete_cases %>%
  mutate(CSAI2R_sum = rowSums(select(., all_of(CSAI2R_items_names)), na.rm = TRUE),
         BE_sum = rowSums(select(., all_of(BE_items_names)), na.rm = TRUE))

#plot of sum scores (raw data, not ranked)
ggplot(data_complete_cases_CSAI2R_BE, aes(x = CSAI2R_sum, y = BE_sum)) +
  geom_point() +
  labs(
    title = "Scatterplot of Performance rating vs BE Sum Scores",
    x = "CSAI-2R Sum Score",
    y = "BE Sum Score"
  ) +
  theme_minimal()

Code
#compute Spearman's rank correlation
correlation_CSAI2R_BE <- correlation(data_complete_cases_CSAI2R_BE[, c("CSAI2R_sum", "BE_sum")], method = "spearman")

print(correlation_CSAI2R_BE)
# Correlation Matrix (spearman-method)

Parameter1 | Parameter2 |   rho |         95% CI |        S |         p
-----------------------------------------------------------------------
CSAI2R_sum |     BE_sum | -0.36 | [-0.52, -0.18] | 3.36e+05 | < .001***

p-value adjustment method: Holm (1979)
Observations: 114

Some participants did not complete BEA and CSAI-2R during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#sensitivity analysis correlation BEA and CSAI-2R
CSAI2R_items_less1day <- data_complete_cases_less_than_1day[, c(
  "CSAI-2R.1",
  "CSAI-2R.2",
  "CSAI-2R.3",
  "CSAI-2R.4",
  "CSAI-2R.5",
  "CSAI-2R.6",
  "CSAI-2R.7",
  "CSAI-2R.8",
  "CSAI-2R.9",
  "CSAI-2R.10",
  "CSAI-2R.11",
  "CSAI-2R.12",
  "CSAI-2R.13",
  "CSAI-2R.14",
  "CSAI-2R.15",
  "CSAI-2R.16",
  "CSAI-2R.17"
)]

BE_items_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

BE_items_less1day_names <- names(BE_items_less1day)
CSAI2R_items_less1day_names <- names(CSAI2R_items_less1day)

#computing sum scores for each scale
data_complete_cases_CSAI2R_BE_less1day <- data_complete_cases_less_than_1day %>%
  mutate(
    CSAI2R_sum = rowSums(select(., all_of(
      CSAI2R_items_less1day_names
    )), na.rm = TRUE),
    BE_sum = rowSums(select(., all_of(
      BE_items_less1day_names
    )), na.rm = TRUE)
  )

#compute Spearman's rank correlation
correlation_CSAI2R_BE_less1day <- correlation(data_complete_cases_CSAI2R_BE_less1day[, c("CSAI2R_sum", "BE_sum")], method = "spearman")

print(correlation_CSAI2R_BE_less1day)
# Correlation Matrix (spearman-method)

Parameter1 | Parameter2 |   rho |         95% CI |        S |         p
-----------------------------------------------------------------------
CSAI2R_sum |     BE_sum | -0.44 | [-0.59, -0.26] | 2.33e+05 | < .001***

p-value adjustment method: Holm (1979)
Observations: 99

Result from sensitivity analysis correlation is somewhat close to the original correlation, at least interpretation wise.

8 Subjective performance rating and BEA validity analysis

Validity between the subjective performance rating and BEA is investigated by correlation analysis. The subjective performance rating has not undergone Rasch analysis (since it is only one item), thus participants sum scores will be correlated. Since sum scores are ordinal data, Spearman’s rank correlation is applied.

Code
#complete cases data frame for validity analyses

#removing missing values in relevant items for complete cases analysis
#only BE items and one item from the other scales in the bulls eye study
#is needed (in this case the SubjektivPrestation item) 
#this is because there are no missing values for the other scales (beside some BE items) 
#after the intervention study (Hockeystudien) participants have been removed
#therefore this data frame will work for all validity analyses 

items_validity_analyses <- c(
  "BE_Competition_recoded",
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded",
  "SubjektivPrestation"
)

data_complete_cases <- bulls_eye_merged_superfinal_data[complete.cases(bulls_eye_merged_superfinal_data[, items_validity_analyses]), ]

#156 participants to begin with
#34 participants from the intervention study (Hockeystudien) were removed (did not complete any scale beside Bulls Eye)
#8 participants (not from the intervention study) with missing values in Bulls Eye items were removed
#results in a complete cases data frame with 114 participants 

complete_cases_dims <- dim_names(data_complete_cases, "Användarnamn") 
count <- dplyr::count
count(complete_cases_dims)
# A tibble: 1 × 1
      n
  <int>
1   114
Code
#item preparation and plot of sum scores
performance_rating <- data_complete_cases[, c("SubjektivPrestation")]

BE_items <- data_complete_cases[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

BE_items_names <- names(BE_items)
performance_rating_names <- names(performance_rating)

#computing sum scores 
data_complete_cases_performance_rating_BE <- data_complete_cases %>%
  mutate(
    performance_rating_sum = rowSums(select(., all_of(
      performance_rating_names
    )), na.rm = TRUE),
    BE_sum = rowSums(select(., all_of(BE_items_names)), na.rm = TRUE)
  )

#plot of sum scores (raw data, not ranked)
ggplot(data_complete_cases_performance_rating_BE, aes(x = performance_rating_sum, y = BE_sum)) +
  geom_point() +
  labs(
    title = "Scatterplot of Performance rating vs BE Sum Scores",
    x = "Performance rating",
    y = "BE Sum Score"
  ) +
  theme_minimal()

Code
#compute Spearman's rank correlation
correlation_performance_rating_BE <- correlation(data_complete_cases_performance_rating_BE[, c("performance_rating_sum", "BE_sum")], method = "spearman")

print(correlation_performance_rating_BE)
# Correlation Matrix (spearman-method)

Parameter1             | Parameter2 |  rho |        95% CI |        S |     p
-----------------------------------------------------------------------------
performance_rating_sum |     BE_sum | 0.17 | [-0.02, 0.35] | 2.04e+05 | 0.063

p-value adjustment method: Holm (1979)
Observations: 114

Some participants did not complete BEA and the subjective performance rating during the same day as intended. Therefore, a sensitivity analysis is conducted.

Code
#creating data frame with only participants that completed Bulls Eye and the other scales during the same day
#since 15 participants did not, 99 participants will remain in the sensitivity analysis
data_complete_cases_less_than_1day <- subset(data_complete_cases, time_latency_survey1_bulls_eye1_days < 1)
Code
#sensitivity correlation analysis BEA and subjective performance rating
performance_rating_less1day <- data_complete_cases_less_than_1day[, c("SubjektivPrestation")]

BE_items_less1day <- data_complete_cases_less_than_1day[, c(
  "BE_Training_recoded",
  "BE_PreparationRecovery_recoded",
  "BE_Life_recoded",
  "BE_Obstacles_Sport_recoded",
  "BE_Obstacles_Life_recoded"
)]

BE_items_less1day_names <- names(BE_items_less1day)
performance_rating_less1day_names <- names(performance_rating_less1day)

#computing sum scores
data_complete_cases_performance_rating_BE_less1day <- data_complete_cases_less_than_1day %>%
  mutate(
    performance_rating_sum = rowSums(select(
      ., all_of(performance_rating_less1day_names)
    ), na.rm = TRUE),
    BE_sum = rowSums(select(., all_of(
      BE_items_less1day_names
    )), na.rm = TRUE)
  )

#compute Spearman's rank correlation
correlation_performance_rating_BE_less1day <- correlation(data_complete_cases_performance_rating_BE_less1day[, c("performance_rating_sum", "BE_sum")], method = "spearman")

print(correlation_performance_rating_BE_less1day)
# Correlation Matrix (spearman-method)

Parameter1             | Parameter2 |  rho |        95% CI |        S |     p
-----------------------------------------------------------------------------
performance_rating_sum |     BE_sum | 0.18 | [-0.03, 0.37] | 1.33e+05 | 0.082

p-value adjustment method: Holm (1979)
Observations: 99

Result from sensitivity analysis correlation is close to the original correlation.